conformalVoronoiMeshZones.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) 2013 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 "conformalVoronoiMesh.H"
27 #include "polyModifyFace.H"
28 #include "polyModifyCell.H"
29 #include "syncTools.H"
30 #include "regionSplit.H"
31 #include "surfaceZonesInfo.H"
32 
33 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
34 
35 void Foam::conformalVoronoiMesh::calcNeighbourCellCentres
36 (
37  const polyMesh& mesh,
38  const pointField& cellCentres,
39  pointField& neiCc
40 ) const
41 {
42  label nBoundaryFaces = mesh.nFaces() - mesh.nInternalFaces();
43 
44  if (neiCc.size() != nBoundaryFaces)
45  {
46  FatalErrorIn("conformalVoronoiMesh::calcNeighbourCellCentres(..)")
47  << "nBoundaries:" << nBoundaryFaces
48  << " neiCc:" << neiCc.size()
49  << abort(FatalError);
50  }
51 
52  const polyBoundaryMesh& patches = mesh.boundaryMesh();
53 
54  forAll(patches, patchI)
55  {
56  const polyPatch& pp = patches[patchI];
57 
58  const labelUList& faceCells = pp.faceCells();
59 
60  label bFaceI = pp.start() - mesh.nInternalFaces();
61 
62  if (pp.coupled())
63  {
64  forAll(faceCells, i)
65  {
66  neiCc[bFaceI] = cellCentres[faceCells[i]];
67  bFaceI++;
68  }
69  }
70  }
71 
72  // Swap coupled boundaries. Apply separation to cc since is coordinate.
74 }
75 
76 
77 void Foam::conformalVoronoiMesh::selectSeparatedCoupledFaces
78 (
79  const polyMesh& mesh,
80  boolList& selected
81 ) const
82 {
83  const polyBoundaryMesh& patches = mesh.boundaryMesh();
84 
85  forAll(patches, patchI)
86  {
87  // Check all coupled. Avoid using .coupled() so we also pick up AMI.
88  if (isA<coupledPolyPatch>(patches[patchI]))
89  {
90  const coupledPolyPatch& cpp = refCast<const coupledPolyPatch>
91  (
92  patches[patchI]
93  );
94 
95  if (cpp.separated() || !cpp.parallel())
96  {
97  forAll(cpp, i)
98  {
99  selected[cpp.start()+i] = true;
100  }
101  }
102  }
103  }
104 }
105 
106 
107 void Foam::conformalVoronoiMesh::findCellZoneInsideWalk
108 (
109  const polyMesh& mesh,
110  const labelList& locationSurfaces, // indices of surfaces with inside point
111  const labelList& faceToSurface, // per face index of named surface
112  labelList& cellToSurface
113 ) const
114 {
115  // Analyse regions. Reuse regionsplit
116  boolList blockedFace(mesh.nFaces());
117  selectSeparatedCoupledFaces(mesh, blockedFace);
118 
119  forAll(faceToSurface, faceI)
120  {
121  if (faceToSurface[faceI] == -1)
122  {
123  blockedFace[faceI] = false;
124  }
125  else
126  {
127  blockedFace[faceI] = true;
128  }
129  }
130  // No need to sync since namedSurfaceIndex already is synced
131 
132  // Set region per cell based on walking
133  regionSplit cellRegion(mesh, blockedFace);
134  blockedFace.clear();
135 
136 
137  // Force calculation of face decomposition (used in findCell)
138  (void)mesh.tetBasePtIs();
139 
140  const PtrList<surfaceZonesInfo>& surfZones =
141  geometryToConformTo().surfZones();
142 
143  // For all locationSurface find the cell
144  forAll(locationSurfaces, i)
145  {
146  label surfI = locationSurfaces[i];
147 
148  const Foam::point& insidePoint = surfZones[surfI].zoneInsidePoint();
149 
150  const word& surfName = geometryToConformTo().geometry().names()[surfI];
151 
152  Info<< " For surface " << surfName
153  << " finding inside point " << insidePoint
154  << endl;
155 
156  // Find the region containing the insidePoint
157  label keepRegionI = -1;
158 
159  label cellI = mesh.findCell(insidePoint);
160 
161  if (cellI != -1)
162  {
163  keepRegionI = cellRegion[cellI];
164  }
165  reduce(keepRegionI, maxOp<label>());
166 
167  Info<< " For surface " << surfName
168  << " found point " << insidePoint << " in cell " << cellI
169  << " in global region " << keepRegionI
170  << " out of " << cellRegion.nRegions() << " regions." << endl;
171 
172  if (keepRegionI == -1)
173  {
175  (
176  "conformalVoronoiMesh::findCellZoneInsideWalk"
177  "(const polyMesh&, const labelList&"
178  ", const labelList&, labelList&)"
179  ) << "Point " << insidePoint
180  << " is not inside the mesh." << nl
181  << "Bounding box of the mesh:" << mesh.bounds()
182  << exit(FatalError);
183  }
184 
185  // Set all cells with this region
186  forAll(cellRegion, cellI)
187  {
188  if (cellRegion[cellI] == keepRegionI)
189  {
190  if (cellToSurface[cellI] == -2)
191  {
192  cellToSurface[cellI] = surfI;
193  }
194  else if (cellToSurface[cellI] != surfI)
195  {
196  WarningIn
197  (
198  "conformalVoronoiMesh::findCellZoneInsideWalk"
199  "(const labelList&, const labelList&"
200  ", const labelList&, const labelList&)"
201  ) << "Cell " << cellI
202  << " at " << mesh.cellCentres()[cellI]
203  << " is inside surface " << surfName
204  << " but already marked as being in zone "
205  << cellToSurface[cellI] << endl
206  << "This can happen if your surfaces are not"
207  << " (sufficiently) closed."
208  << endl;
209  }
210  }
211  }
212  }
213 }
214 
215 
216 Foam::labelList Foam::conformalVoronoiMesh::calcCellZones
217 (
218  const pointField& cellCentres
219 ) const
220 {
221  labelList cellToSurface(cellCentres.size(), -1);
222 
223  const PtrList<surfaceZonesInfo>& surfZones =
224  geometryToConformTo().surfZones();
225 
226  // Get list of closed surfaces
227  labelList closedNamedSurfaces
228  (
230  (
231  surfZones,
232  geometryToConformTo().geometry(),
233  geometryToConformTo().surfaces()
234  )
235  );
236 
237  forAll(closedNamedSurfaces, i)
238  {
239  label surfI = closedNamedSurfaces[i];
240 
241  const searchableSurface& surface =
242  allGeometry()[geometryToConformTo().surfaces()[surfI]];
243 
244  const surfaceZonesInfo::areaSelectionAlgo selectionMethod =
245  surfZones[surfI].zoneInside();
246 
247  if
248  (
249  selectionMethod != surfaceZonesInfo::INSIDE
250  && selectionMethod != surfaceZonesInfo::OUTSIDE
251  && selectionMethod != surfaceZonesInfo::INSIDEPOINT
252  )
253  {
254  FatalErrorIn("conformalVoronoiMesh::calcCellZones(..)")
255  << "Trying to use surface "
256  << surface.name()
257  << " which has non-geometric inside selection method "
258  << surfaceZonesInfo::areaSelectionAlgoNames[selectionMethod]
259  << exit(FatalError);
260  }
261 
262  if (surface.hasVolumeType())
263  {
264  List<volumeType> volType;
265  surface.getVolumeType(cellCentres, volType);
266 
267  bool selectInside = true;
268  if (selectionMethod == surfaceZonesInfo::INSIDEPOINT)
269  {
270  List<volumeType> volTypeInsidePoint;
271  surface.getVolumeType
272  (
273  pointField(1, surfZones[surfI].zoneInsidePoint()),
274  volTypeInsidePoint
275  );
276 
277  if (volTypeInsidePoint[0] == volumeType::OUTSIDE)
278  {
279  selectInside = false;
280  }
281  }
282  else if (selectionMethod == surfaceZonesInfo::OUTSIDE)
283  {
284  selectInside = false;
285  }
286 
287  forAll(volType, pointI)
288  {
289  if (cellToSurface[pointI] == -1)
290  {
291  if
292  (
293  (
294  volType[pointI] == volumeType::INSIDE
295  && selectInside
296  )
297  || (
298  volType[pointI] == volumeType::OUTSIDE
299  && !selectInside
300  )
301  )
302  {
303  cellToSurface[pointI] = surfI;
304  }
305  }
306  }
307  }
308  }
309 
310  return cellToSurface;
311 }
312 
313 
314 void Foam::conformalVoronoiMesh::calcFaceZones
315 (
316  const polyMesh& mesh,
317  const pointField& cellCentres,
318  const labelList& cellToSurface,
319  labelList& faceToSurface,
320  boolList& flipMap
321 ) const
322 {
323  faceToSurface.setSize(mesh.nFaces(), -1);
324  flipMap.setSize(mesh.nFaces(), false);
325 
326  const faceList& faces = mesh.faces();
327  const labelList& faceOwner = mesh.faceOwner();
328  const labelList& faceNeighbour = mesh.faceNeighbour();
329 
330  labelList neiFaceOwner(mesh.nFaces() - mesh.nInternalFaces(), -1);
331 
332  const polyBoundaryMesh& patches = mesh.boundaryMesh();
333 
334  forAll(patches, patchI)
335  {
336  const polyPatch& pp = patches[patchI];
337 
338  const labelUList& faceCells = pp.faceCells();
339 
340  label bFaceI = pp.start() - mesh.nInternalFaces();
341 
342  if (pp.coupled())
343  {
344  forAll(faceCells, i)
345  {
346  neiFaceOwner[bFaceI] = cellToSurface[faceCells[i]];
347  bFaceI++;
348  }
349  }
350  }
351 
352  syncTools::swapBoundaryFaceList(mesh, neiFaceOwner);
353 
354  forAll(faces, faceI)
355  {
356  const label ownerSurfaceI = cellToSurface[faceOwner[faceI]];
357 
358  if (faceToSurface[faceI] >= 0)
359  {
360  continue;
361  }
362 
363  if (mesh.isInternalFace(faceI))
364  {
365  const label neiSurfaceI = cellToSurface[faceNeighbour[faceI]];
366 
367  if
368  (
369  (ownerSurfaceI >= 0 || neiSurfaceI >= 0)
370  && ownerSurfaceI != neiSurfaceI
371  )
372  {
373  flipMap[faceI] =
374  (
375  ownerSurfaceI == max(ownerSurfaceI, neiSurfaceI)
376  ? false
377  : true
378  );
379 
380  faceToSurface[faceI] = max(ownerSurfaceI, neiSurfaceI);
381  }
382  }
383  else
384  {
385  label patchID = mesh.boundaryMesh().whichPatch(faceI);
386 
387  if (mesh.boundaryMesh()[patchID].coupled())
388  {
389  const label neiSurfaceI =
390  neiFaceOwner[faceI - mesh.nInternalFaces()];
391 
392  if
393  (
394  (ownerSurfaceI >= 0 || neiSurfaceI >= 0)
395  && ownerSurfaceI != neiSurfaceI
396  )
397  {
398  flipMap[faceI] =
399  (
400  ownerSurfaceI == max(ownerSurfaceI, neiSurfaceI)
401  ? false
402  : true
403  );
404 
405  faceToSurface[faceI] = max(ownerSurfaceI, neiSurfaceI);
406  }
407  }
408  else
409  {
410  if (ownerSurfaceI >= 0)
411  {
412  faceToSurface[faceI] = ownerSurfaceI;
413  }
414  }
415  }
416  }
417 
418 
419  const PtrList<surfaceZonesInfo>& surfZones =
420  geometryToConformTo().surfZones();
421 
422  labelList unclosedSurfaces
423  (
425  (
426  surfZones,
427  geometryToConformTo().geometry(),
428  geometryToConformTo().surfaces()
429  )
430  );
431 
432  pointField neiCc(mesh.nFaces() - mesh.nInternalFaces());
433  calcNeighbourCellCentres
434  (
435  mesh,
436  cellCentres,
437  neiCc
438  );
439 
440  // Use intersection of cellCentre connections
441  forAll(faces, faceI)
442  {
443  if (faceToSurface[faceI] >= 0)
444  {
445  continue;
446  }
447 
448  label patchID = mesh.boundaryMesh().whichPatch(faceI);
449 
450  const label own = faceOwner[faceI];
451 
452  List<pointIndexHit> surfHit;
453  labelList hitSurface;
454 
455  if (mesh.isInternalFace(faceI))
456  {
457  const label nei = faceNeighbour[faceI];
458 
459  geometryToConformTo().findSurfaceAllIntersections
460  (
461  cellCentres[own],
462  cellCentres[nei],
463  surfHit,
464  hitSurface
465  );
466  }
467  else if (patchID != -1 && mesh.boundaryMesh()[patchID].coupled())
468  {
469  geometryToConformTo().findSurfaceAllIntersections
470  (
471  cellCentres[own],
472  neiCc[faceI - mesh.nInternalFaces()],
473  surfHit,
474  hitSurface
475  );
476  }
477 
478  // If there are multiple intersections then do not add to
479  // a faceZone
480  if (surfHit.size() == 1 && surfHit[0].hit())
481  {
482  if (findIndex(unclosedSurfaces, hitSurface[0]) != -1)
483  {
484  vectorField norm;
485  geometryToConformTo().getNormal
486  (
487  hitSurface[0],
488  List<pointIndexHit>(1, surfHit[0]),
489  norm
490  );
491 
492  vector fN = faces[faceI].normal(mesh.points());
493  fN /= mag(fN) + SMALL;
494 
495  if ((norm[0] & fN) < 0)
496  {
497  flipMap[faceI] = true;
498  }
499  else
500  {
501  flipMap[faceI] = false;
502  }
503 
504  faceToSurface[faceI] = hitSurface[0];
505  }
506  }
507  }
508 
509 
510 // labelList neiCellSurface(mesh.nFaces()-mesh.nInternalFaces());
511 //
512 // forAll(patches, patchI)
513 // {
514 // const polyPatch& pp = patches[patchI];
515 //
516 // if (pp.coupled())
517 // {
518 // forAll(pp, i)
519 // {
520 // label faceI = pp.start()+i;
521 // label ownSurface = cellToSurface[faceOwner[faceI]];
522 // neiCellSurface[faceI - mesh.nInternalFaces()] = ownSurface;
523 // }
524 // }
525 // }
526 // syncTools::swapBoundaryFaceList(mesh, neiCellSurface);
527 //
528 // forAll(patches, patchI)
529 // {
530 // const polyPatch& pp = patches[patchI];
531 //
532 // if (pp.coupled())
533 // {
534 // forAll(pp, i)
535 // {
536 // label faceI = pp.start()+i;
537 // label ownSurface = cellToSurface[faceOwner[faceI]];
538 // label neiSurface =
539 // neiCellSurface[faceI-mesh.nInternalFaces()];
540 //
541 // if (faceToSurface[faceI] == -1 && (ownSurface != neiSurface))
542 // {
543 // // Give face the max cell zone
544 // faceToSurface[faceI] = max(ownSurface, neiSurface);
545 // }
546 // }
547 // }
548 // }
549 
550  // Sync
551  syncTools::syncFaceList(mesh, faceToSurface, maxEqOp<label>());
552 }
553 
554 
555 void Foam::conformalVoronoiMesh::addZones
556 (
557  polyMesh& mesh,
558  const pointField& cellCentres
559 ) const
560 {
561  Info<< " Adding zones to mesh" << endl;
562 
563  const PtrList<surfaceZonesInfo>& surfZones =
564  geometryToConformTo().surfZones();
565 
566  labelList cellToSurface(calcCellZones(cellCentres));
567 
568  labelList faceToSurface;
569  boolList flipMap;
570 
571  calcFaceZones
572  (
573  mesh,
574  cellCentres,
575  cellToSurface,
576  faceToSurface,
577  flipMap
578  );
579 
580  labelList insidePointNamedSurfaces
581  (
583  );
584 
585  findCellZoneInsideWalk
586  (
587  mesh,
588  insidePointNamedSurfaces,
589  faceToSurface,
590  cellToSurface
591  );
592 
593  labelList namedSurfaces(surfaceZonesInfo::getNamedSurfaces(surfZones));
594 
595  forAll(namedSurfaces, i)
596  {
597  label surfI = namedSurfaces[i];
598 
599  Info<< incrIndent << indent << "Surface : "
600  << geometryToConformTo().geometry().names()[surfI] << nl
601  << indent << " faceZone : "
602  << surfZones[surfI].faceZoneName() << nl
603  << indent << " cellZone : "
604  << surfZones[surfI].cellZoneName()
605  << decrIndent << endl;
606  }
607 
608  // Add zones to mesh
609  labelList surfaceToFaceZone =
611  (
612  surfZones,
613  namedSurfaces,
614  mesh
615  );
616 
617  labelList surfaceToCellZone =
619  (
620  surfZones,
621  namedSurfaces,
622  mesh
623  );
624 
625  // Topochange container
626  polyTopoChange meshMod(mesh);
627 
628  forAll(cellToSurface, cellI)
629  {
630  label surfaceI = cellToSurface[cellI];
631 
632  if (surfaceI >= 0)
633  {
634  label zoneI = surfaceToCellZone[surfaceI];
635 
636  if (zoneI >= 0)
637  {
638  meshMod.setAction
639  (
640  polyModifyCell
641  (
642  cellI,
643  false, // removeFromZone
644  zoneI
645  )
646  );
647  }
648  }
649  }
650 
651  const labelList& faceOwner = mesh.faceOwner();
652  const labelList& faceNeighbour = mesh.faceNeighbour();
653 
654  forAll(faceToSurface, faceI)
655  {
656  label surfaceI = faceToSurface[faceI];
657 
658  if (surfaceI < 0)
659  {
660  continue;
661  }
662 
663  label patchID = mesh.boundaryMesh().whichPatch(faceI);
664 
665  if (mesh.isInternalFace(faceI))
666  {
667  label own = faceOwner[faceI];
668  label nei = faceNeighbour[faceI];
669 
670  meshMod.setAction
671  (
672  polyModifyFace
673  (
674  mesh.faces()[faceI], // modified face
675  faceI, // label of face
676  own, // owner
677  nei, // neighbour
678  false, // face flip
679  -1, // patch for face
680  false, // remove from zone
681  surfaceToFaceZone[surfaceI], // zone for face
682  flipMap[faceI] // face flip in zone
683  )
684  );
685  }
686  else if (patchID != -1 && mesh.boundaryMesh()[patchID].coupled())
687  {
688  label own = faceOwner[faceI];
689 
690  meshMod.setAction
691  (
692  polyModifyFace
693  (
694  mesh.faces()[faceI], // modified face
695  faceI, // label of face
696  own, // owner
697  -1, // neighbour
698  false, // face flip
699  patchID, // patch for face
700  false, // remove from zone
701  surfaceToFaceZone[surfaceI], // zone for face
702  flipMap[faceI] // face flip in zone
703  )
704  );
705  }
706  }
707 
708  // Change the mesh (no inflation, parallel sync)
709  autoPtr<mapPolyMesh> map = meshMod.changeMesh(mesh, false, true);
710 }
711 
712 
713 // ************************************************************************* //
static void syncFaceList(const polyMesh &mesh, UList< T > &l, const CombineOp &cop)
Synchronize values on all mesh faces.
Definition: syncTools.H:380
Field< vector > vectorField
Specialisation of Field<T> for vector.
static labelList getInsidePointNamedSurfaces(const PtrList< surfaceZonesInfo > &surfList)
Get indices of surfaces with a cellZone that have 'insidePoint'.
List< label > labelList
A List of labels.
Definition: labelList.H:56
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:248
void reduce(const List< UPstream::commsStruct > &comms, T &Value, const BinaryOp &bop, const int tag, const label comm)
error FatalError
static const NamedEnum< areaSelectionAlgo, 4 > areaSelectionAlgoNames
dimensioned< scalar > mag(const dimensioned< Type > &)
Ostream & indent(Ostream &os)
Indent stream.
Definition: Ostream.H:218
Ostream & incrIndent(Ostream &os)
Increment the indent level.
Definition: Ostream.H:225
UList< label > labelUList
Definition: UList.H:62
Ostream & decrIndent(Ostream &os)
Decrement the indent level.
Definition: Ostream.H:232
#define WarningIn(fn)
static void swapBoundaryFaceList(const polyMesh &mesh, UList< T > &l)
Swap coupled boundary face values.
Definition: syncTools.H:429
static labelList addFaceZonesToMesh(const PtrList< surfaceZonesInfo > &surfList, const labelList &namedSurfaces, polyMesh &mesh)
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
label findIndex(const ListType &, typename ListType::const_reference, const label start=0)
Find first occurence of given element and return index,.
static labelList getUnclosedNamedSurfaces(const PtrList< surfaceZonesInfo > &surfList, const searchableSurfaces &allGeometry, const labelList &surfaces)
Get indices of surfaces with a cellZone that are unclosed.
errorManip< error > abort(error &err)
Definition: errorManip.H:131
messageStream Info
static labelList getAllClosedNamedSurfaces(const PtrList< surfaceZonesInfo > &surfList, const searchableSurfaces &allGeometry, const labelList &surfaces)
Get indices of surfaces with a cellZone that are closed.
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:124
areaSelectionAlgo
Types of selection of area.
#define forAll(list, i)
Definition: UList.H:420
#define FatalErrorIn(fn)
Definition: error.H:317
static const char nl
Definition: Ostream.H:257
vectorField pointField
pointField is a vectorField.
Definition: pointFieldFwd.H:42
static labelList getNamedSurfaces(const PtrList< surfaceZonesInfo > &surfList)
Get indices of named surfaces (surfaces with faceZoneName)
static void swapBoundaryFacePositions(const polyMesh &mesh, UList< point > &l)
Swap coupled positions.
Definition: syncTools.H:445
List< face > faceList
Definition: faceListFwd.H:43
static labelList addCellZonesToMesh(const PtrList< surfaceZonesInfo > &surfList, const labelList &namedSurfaces, polyMesh &mesh)
Vector< scalar > vector
A scalar version of the templated Vector.
Definition: vector.H:49
List< bool > boolList
Bool container classes.
Definition: boolList.H:50