conformalVoronoiMeshZones.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) 2013-2018 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  {
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 =
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  << "Point " << insidePoint
176  << " is not inside the mesh." << nl
177  << "Bounding box of the mesh:" << mesh.bounds()
178  << exit(FatalError);
179  }
180 
181  // Set all cells with this region
182  forAll(cellRegion, celli)
183  {
184  if (cellRegion[celli] == keepRegionI)
185  {
186  if (cellToSurface[celli] == -2)
187  {
188  cellToSurface[celli] = surfI;
189  }
190  else if (cellToSurface[celli] != surfI)
191  {
193  << "Cell " << celli
194  << " at " << mesh.cellCentres()[celli]
195  << " is inside surface " << surfName
196  << " but already marked as being in zone "
197  << cellToSurface[celli] << endl
198  << "This can happen if your surfaces are not"
199  << " (sufficiently) closed."
200  << endl;
201  }
202  }
203  }
204  }
205 }
206 
207 
208 Foam::labelList Foam::conformalVoronoiMesh::calcCellZones
209 (
210  const pointField& cellCentres
211 ) const
212 {
213  labelList cellToSurface(cellCentres.size(), label(-1));
214 
215  const PtrList<surfaceZonesInfo>& surfZones =
217 
218  // Get list of closed surfaces
219  labelList closedNamedSurfaces
220  (
222  (
223  surfZones,
224  geometryToConformTo().geometry(),
225  geometryToConformTo().surfaces()
226  )
227  );
228 
229  forAll(closedNamedSurfaces, i)
230  {
231  label surfI = closedNamedSurfaces[i];
232 
233  const searchableSurface& surface =
235 
236  const surfaceZonesInfo::areaSelectionAlgo selectionMethod =
237  surfZones[surfI].zoneInside();
238 
239  if
240  (
241  selectionMethod != surfaceZonesInfo::INSIDE
242  && selectionMethod != surfaceZonesInfo::OUTSIDE
243  && selectionMethod != surfaceZonesInfo::INSIDEPOINT
244  )
245  {
247  << "Trying to use surface "
248  << surface.name()
249  << " which has non-geometric inside selection method "
250  << surfaceZonesInfo::areaSelectionAlgoNames[selectionMethod]
251  << exit(FatalError);
252  }
253 
254  if (surface.hasVolumeType())
255  {
256  List<volumeType> volType;
257  surface.getVolumeType(cellCentres, volType);
258 
259  bool selectInside = true;
260  if (selectionMethod == surfaceZonesInfo::INSIDEPOINT)
261  {
262  List<volumeType> volTypeInsidePoint;
263  surface.getVolumeType
264  (
265  pointField(1, surfZones[surfI].zoneInsidePoint()),
266  volTypeInsidePoint
267  );
268 
269  if (volTypeInsidePoint[0] == volumeType::OUTSIDE)
270  {
271  selectInside = false;
272  }
273  }
274  else if (selectionMethod == surfaceZonesInfo::OUTSIDE)
275  {
276  selectInside = false;
277  }
278 
279  forAll(volType, pointi)
280  {
281  if (cellToSurface[pointi] == -1)
282  {
283  if
284  (
285  (
286  volType[pointi] == volumeType::INSIDE
287  && selectInside
288  )
289  || (
290  volType[pointi] == volumeType::OUTSIDE
291  && !selectInside
292  )
293  )
294  {
295  cellToSurface[pointi] = surfI;
296  }
297  }
298  }
299  }
300  }
301 
302  return cellToSurface;
303 }
304 
305 
306 void Foam::conformalVoronoiMesh::calcFaceZones
307 (
308  const polyMesh& mesh,
309  const pointField& cellCentres,
310  const labelList& cellToSurface,
311  labelList& faceToSurface,
312  boolList& flipMap
313 ) const
314 {
315  faceToSurface.setSize(mesh.nFaces(), -1);
316  flipMap.setSize(mesh.nFaces(), false);
317 
318  const faceList& faces = mesh.faces();
319  const labelList& faceOwner = mesh.faceOwner();
320  const labelList& faceNeighbour = mesh.faceNeighbour();
321 
322  labelList neiFaceOwner(mesh.nFaces() - mesh.nInternalFaces(), label(-1));
323 
324  const polyBoundaryMesh& patches = mesh.boundaryMesh();
325 
326  forAll(patches, patchi)
327  {
328  const polyPatch& pp = patches[patchi];
329 
330  const labelUList& faceCells = pp.faceCells();
331 
332  label bFacei = pp.start() - mesh.nInternalFaces();
333 
334  if (pp.coupled())
335  {
336  forAll(faceCells, i)
337  {
338  neiFaceOwner[bFacei] = cellToSurface[faceCells[i]];
339  bFacei++;
340  }
341  }
342  }
343 
344  syncTools::swapBoundaryFaceList(mesh, neiFaceOwner);
345 
346  forAll(faces, facei)
347  {
348  const label ownerSurfacei = cellToSurface[faceOwner[facei]];
349 
350  if (faceToSurface[facei] >= 0)
351  {
352  continue;
353  }
354 
355  if (mesh.isInternalFace(facei))
356  {
357  const label neiSurfacei = cellToSurface[faceNeighbour[facei]];
358 
359  if
360  (
361  (ownerSurfacei >= 0 || neiSurfacei >= 0)
362  && ownerSurfacei != neiSurfacei
363  )
364  {
365  flipMap[facei] =
366  (
367  ownerSurfacei == max(ownerSurfacei, neiSurfacei)
368  ? false
369  : true
370  );
371 
372  faceToSurface[facei] = max(ownerSurfacei, neiSurfacei);
373  }
374  }
375  else
376  {
377  label patchID = mesh.boundaryMesh().whichPatch(facei);
378 
379  if (mesh.boundaryMesh()[patchID].coupled())
380  {
381  const label neiSurfacei =
382  neiFaceOwner[facei - mesh.nInternalFaces()];
383 
384  if
385  (
386  (ownerSurfacei >= 0 || neiSurfacei >= 0)
387  && ownerSurfacei != neiSurfacei
388  )
389  {
390  flipMap[facei] =
391  (
392  ownerSurfacei == max(ownerSurfacei, neiSurfacei)
393  ? false
394  : true
395  );
396 
397  faceToSurface[facei] = max(ownerSurfacei, neiSurfacei);
398  }
399  }
400  else
401  {
402  if (ownerSurfacei >= 0)
403  {
404  faceToSurface[facei] = ownerSurfacei;
405  }
406  }
407  }
408  }
409 
410 
411  const PtrList<surfaceZonesInfo>& surfZones =
413 
414  labelList unclosedSurfaces
415  (
417  (
418  surfZones,
419  geometryToConformTo().geometry(),
420  geometryToConformTo().surfaces()
421  )
422  );
423 
424  pointField neiCc(mesh.nFaces() - mesh.nInternalFaces());
425  calcNeighbourCellCentres
426  (
427  mesh,
428  cellCentres,
429  neiCc
430  );
431 
432  // Use intersection of cellCentre connections
433  forAll(faces, facei)
434  {
435  if (faceToSurface[facei] >= 0)
436  {
437  continue;
438  }
439 
440  label patchID = mesh.boundaryMesh().whichPatch(facei);
441 
442  const label own = faceOwner[facei];
443 
444  List<pointIndexHit> surfHit;
445  labelList hitSurface;
446 
447  if (mesh.isInternalFace(facei))
448  {
449  const label nei = faceNeighbour[facei];
450 
452  (
453  cellCentres[own],
454  cellCentres[nei],
455  surfHit,
456  hitSurface
457  );
458  }
459  else if (patchID != -1 && mesh.boundaryMesh()[patchID].coupled())
460  {
462  (
463  cellCentres[own],
464  neiCc[facei - mesh.nInternalFaces()],
465  surfHit,
466  hitSurface
467  );
468  }
469 
470  // If there are multiple intersections then do not add to
471  // a faceZone
472  if (surfHit.size() == 1 && surfHit[0].hit())
473  {
474  if (findIndex(unclosedSurfaces, hitSurface[0]) != -1)
475  {
476  vectorField norm;
478  (
479  hitSurface[0],
480  List<pointIndexHit>(1, surfHit[0]),
481  norm
482  );
483 
484  vector fN = faces[facei].area(mesh.points());
485  fN /= mag(fN) + small;
486 
487  if ((norm[0] & fN) < 0)
488  {
489  flipMap[facei] = true;
490  }
491  else
492  {
493  flipMap[facei] = false;
494  }
495 
496  faceToSurface[facei] = hitSurface[0];
497  }
498  }
499  }
500 
501 
502 // labelList neiCellSurface(mesh.nFaces()-mesh.nInternalFaces());
503 //
504 // forAll(patches, patchi)
505 // {
506 // const polyPatch& pp = patches[patchi];
507 //
508 // if (pp.coupled())
509 // {
510 // forAll(pp, i)
511 // {
512 // label facei = pp.start()+i;
513 // label ownSurface = cellToSurface[faceOwner[facei]];
514 // neiCellSurface[facei - mesh.nInternalFaces()] = ownSurface;
515 // }
516 // }
517 // }
518 // syncTools::swapBoundaryFaceList(mesh, neiCellSurface);
519 //
520 // forAll(patches, patchi)
521 // {
522 // const polyPatch& pp = patches[patchi];
523 //
524 // if (pp.coupled())
525 // {
526 // forAll(pp, i)
527 // {
528 // label facei = pp.start()+i;
529 // label ownSurface = cellToSurface[faceOwner[facei]];
530 // label neiSurface =
531 // neiCellSurface[facei-mesh.nInternalFaces()];
532 //
533 // if (faceToSurface[facei] == -1 && (ownSurface != neiSurface))
534 // {
535 // // Give face the max cell zone
536 // faceToSurface[facei] = max(ownSurface, neiSurface);
537 // }
538 // }
539 // }
540 // }
541 
542  // Sync
543  syncTools::syncFaceList(mesh, faceToSurface, maxEqOp<label>());
544 }
545 
546 
547 void Foam::conformalVoronoiMesh::addZones
548 (
549  polyMesh& mesh,
550  const pointField& cellCentres
551 ) const
552 {
553  Info<< " Adding zones to mesh" << endl;
554 
555  const PtrList<surfaceZonesInfo>& surfZones =
557 
558  labelList cellToSurface(calcCellZones(cellCentres));
559 
560  labelList faceToSurface;
561  boolList flipMap;
562 
563  calcFaceZones
564  (
565  mesh,
566  cellCentres,
567  cellToSurface,
568  faceToSurface,
569  flipMap
570  );
571 
572  labelList insidePointNamedSurfaces
573  (
575  );
576 
577  findCellZoneInsideWalk
578  (
579  mesh,
580  insidePointNamedSurfaces,
581  faceToSurface,
582  cellToSurface
583  );
584 
585  labelList namedSurfaces(surfaceZonesInfo::getNamedSurfaces(surfZones));
586 
587  forAll(namedSurfaces, i)
588  {
589  label surfI = namedSurfaces[i];
590 
591  Info<< incrIndent << indent << "Surface : "
592  << geometryToConformTo().geometry().names()[surfI] << nl
593  << indent << " faceZone : "
594  << surfZones[surfI].faceZoneName() << nl
595  << indent << " cellZone : "
596  << surfZones[surfI].cellZoneName()
597  << decrIndent << endl;
598  }
599 
600  // Add zones to mesh
601  labelList surfaceToFaceZone =
603  (
604  surfZones,
605  namedSurfaces,
606  mesh
607  );
608 
609  labelList surfaceToCellZone =
611  (
612  surfZones,
613  namedSurfaces,
614  mesh
615  );
616 
617  // Topochange container
618  polyTopoChange meshMod(mesh);
619 
620  forAll(cellToSurface, celli)
621  {
622  label surfacei = cellToSurface[celli];
623 
624  if (surfacei >= 0)
625  {
626  label zoneI = surfaceToCellZone[surfacei];
627 
628  if (zoneI >= 0)
629  {
630  meshMod.setAction
631  (
632  polyModifyCell
633  (
634  celli,
635  false, // removeFromZone
636  zoneI
637  )
638  );
639  }
640  }
641  }
642 
643  const labelList& faceOwner = mesh.faceOwner();
644  const labelList& faceNeighbour = mesh.faceNeighbour();
645 
646  forAll(faceToSurface, facei)
647  {
648  label surfacei = faceToSurface[facei];
649 
650  if (surfacei < 0)
651  {
652  continue;
653  }
654 
655  label patchID = mesh.boundaryMesh().whichPatch(facei);
656 
657  if (mesh.isInternalFace(facei))
658  {
659  label own = faceOwner[facei];
660  label nei = faceNeighbour[facei];
661 
662  meshMod.setAction
663  (
664  polyModifyFace
665  (
666  mesh.faces()[facei], // modified face
667  facei, // label of face
668  own, // owner
669  nei, // neighbour
670  false, // face flip
671  -1, // patch for face
672  false, // remove from zone
673  surfaceToFaceZone[surfacei], // zone for face
674  flipMap[facei] // face flip in zone
675  )
676  );
677  }
678  else if (patchID != -1 && mesh.boundaryMesh()[patchID].coupled())
679  {
680  label own = faceOwner[facei];
681 
682  meshMod.setAction
683  (
684  polyModifyFace
685  (
686  mesh.faces()[facei], // modified face
687  facei, // label of face
688  own, // owner
689  -1, // neighbour
690  false, // face flip
691  patchID, // patch for face
692  false, // remove from zone
693  surfaceToFaceZone[surfacei], // zone for face
694  flipMap[facei] // face flip in zone
695  )
696  );
697  }
698  }
699 
700  // Change the mesh (no inflation, parallel sync)
701  autoPtr<mapPolyMesh> map = meshMod.changeMesh(mesh, false, true);
702 }
703 
704 
705 // ************************************************************************* //
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:428
const searchableSurfaces & geometry() const
Return reference to the searchableSurfaces object containing all.
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
Ostream & indent(Ostream &os)
Indent stream.
Definition: Ostream.H:226
const searchableSurfaces & allGeometry() const
Return the allGeometry object.
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:124
error FatalError
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:319
static void swapBoundaryFacePositions(const polyMesh &mesh, UList< point > &l)
Swap coupled positions.
Definition: syncTools.H:446
const labelList & surfaces() const
Return the surface indices.
List< face > faceList
Definition: faceListFwd.H:43
const PtrList< surfaceZonesInfo > & surfZones() const
Return the surfaceZonesInfo.
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:256
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
const conformationSurfaces & geometryToConformTo() const
Return the conformationSurfaces object.
static labelList getUnclosedNamedSurfaces(const PtrList< surfaceZonesInfo > &surfList, const searchableSurfaces &allGeometry, const labelList &surfaces)
Get indices of surfaces with a cellZone that are unclosed.
UList< label > labelUList
Definition: UList.H:64
static labelList getAllClosedNamedSurfaces(const PtrList< surfaceZonesInfo > &surfList, const searchableSurfaces &allGeometry, const labelList &surfaces)
Get indices of surfaces with a cellZone that are closed.
List< bool > boolList
Bool container classes.
Definition: boolList.H:50
vectorField pointField
pointField is a vectorField.
Definition: pointFieldFwd.H:42
void getNormal(const label hitSurface, const List< pointIndexHit > &surfHit, vectorField &normal) const
List< label > labelList
A List of labels.
Definition: labelList.H:56
static labelList getInsidePointNamedSurfaces(const PtrList< surfaceZonesInfo > &surfList)
Get indices of surfaces with a cellZone that have &#39;insidePoint&#39;.
errorManip< error > abort(error &err)
Definition: errorManip.H:131
static const char nl
Definition: Ostream.H:265
Ostream & decrIndent(Ostream &os)
Decrement the indent level.
Definition: Ostream.H:240
label findIndex(const ListType &, typename ListType::const_reference, const label start=0)
Find first occurrence of given element and return index,.
void reduce(const List< UPstream::commsStruct > &comms, T &Value, const BinaryOp &bop, const int tag, const label comm)
void findSurfaceAllIntersections(const point &start, const point &end, List< pointIndexHit > &surfHit, labelList &hitSurface) const
static void syncFaceList(const polyMesh &mesh, UList< T > &l, const CombineOp &cop)
Synchronize values on all mesh faces.
Definition: syncTools.H:381
label patchi
#define WarningInFunction
Report a warning using Foam::Warning.
const wordList & names() const
static labelList addFaceZonesToMesh(const PtrList< surfaceZonesInfo > &surfList, const labelList &namedSurfaces, polyMesh &mesh)
static labelList getNamedSurfaces(const PtrList< surfaceZonesInfo > &surfList)
Get indices of named surfaces (surfaces with faceZoneName)
messageStream Info
static const NamedEnum< areaSelectionAlgo, 4 > areaSelectionAlgoNames
dimensioned< scalar > mag(const dimensioned< Type > &)
Field< vector > vectorField
Specialisation of Field<T> for vector.
static void swapBoundaryFaceList(const polyMesh &mesh, UList< T > &l)
Swap coupled boundary face values.
Definition: syncTools.H:430
Ostream & incrIndent(Ostream &os)
Increment the indent level.
Definition: Ostream.H:233
areaSelectionAlgo
Types of selection of area.