polyMeshUnMergeCyclics.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) 2020-2022 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 "polyMeshUnMergeCyclics.H"
27 #include "cyclicPolyPatch.H"
28 #include "mergedCyclicPolyPatch.H"
29 #include "polyTopoChange.H"
30 #include "unitConversion.H"
31 
32 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
33 
34 namespace Foam
35 {
36 
38 (
39  const primitivePatch& pp,
40  const scalar includedAngle
41 )
42 {
43  const scalar minCos = Foam::cos(degToRad(180.0 - includedAngle));
44 
45  const labelListList& edgeFaces = pp.edgeFaces();
46  const labelListList& faceEdges = pp.faceEdges();
47  const vectorField& faceNormals = pp.faceNormals();
48 
49  // Working arrays
50  labelList facePrev(pp.size(), -1);
51  labelList faceEdge0(pp.size(), -1);
52  labelList faceEdgeFace0(pp.size(), -1);
53 
54  // The resulting set of zones
55  labelList faceZones(pp.size(), -1);
56 
57  // The current zone being populated
58  label zonei = 0;
59 
60  // The starting face of the current zone
61  label facei0 = 0;
62 
63  // While there are faces to consider
64  while (facei0 < faceZones.size())
65  {
66  label facei = facei0;
67  label faceEdgei = 0;
68  label faceEdgeFacei = 0;
69 
70  faceZones[facei] = zonei;
71 
72  do
73  {
74  // Find the connected edge and face
75  const label edgei = faceEdges[facei][faceEdgei];
76  const label facej = edgeFaces[edgei][faceEdgeFacei];
77 
78  // Find the corresponding position within the other face
79  const label faceEdgej = findIndex(faceEdges[facej], edgei);
80  const label faceEdgeFacej = findIndex(edgeFaces[edgei], facej);
81 
82  if
83  (
84  faceZones[facej] == zonei
85  || (faceNormals[facei] & faceNormals[facej]) < minCos
86  || face::compare(pp[facei], pp[facej]) == -1
87  || edge::compare
88  (
89  pp[facei].faceEdge(faceEdgei),
90  pp[facej].faceEdge(faceEdgej)
91  ) != -1
92  )
93  {
94  // Move to the next face for this edge
95  faceEdgeFacei = edgeFaces[edgei].fcIndex(faceEdgeFacei);
96 
97  // If wrapped around, move to the next edge on the face
98  if (faceEdgeFacei == 0)
99  {
100  faceEdgei = faceEdges[facei].fcIndex(faceEdgei);
101  }
102 
103  // If back at the start then backtrack to the previous face
104  if
105  (
106  faceEdgei == faceEdge0[facei]
107  && faceEdgeFacei == faceEdgeFace0[facei]
108  )
109  {
110  facei = facePrev[facei];
111  faceEdgei = faceEdge0[facei];
112  faceEdgeFacei = faceEdgeFace0[facei];
113  }
114  }
115  else
116  {
117  // Move into the next face
118  faceZones[facej] = zonei;
119  facePrev[facej] = facei;
120  facei = facej;
121 
122  // Set the starting position within this face
123  faceEdgei = faceEdge0[facej] = faceEdgej;
124  faceEdgeFacei = faceEdgeFace0[facej] = faceEdgeFacej;
125  }
126  }
127  while (facei != facei0 || faceEdgei != 0 || faceEdgeFacei != 0);
128 
129  // Get the next zone and starting face
130  ++ zonei;
131  while (facei0 < faceZones.size() && faceZones[facei0] != -1)
132  {
133  ++ facei0;
134  }
135  }
136 
137  return faceZones;
138 }
139 
140 
142 (
143  const primitivePatch& pp,
144  const scalar includedAngle
145 )
146 {
147  const labelList faceZones(primitivePatchGetZones(pp, includedAngle));
148 
149  // If the number of zones is not two, then the split has failed
150  // !!! To improve the robustness of this in the presence of non-planar
151  // cyclics we would have to find pairs of zones that all share the same
152  // transformation. We would still fail if the number of zones were odd.
153  if (findMin(faceZones) != 0 && findMax(faceZones) != 1)
154  {
156  << "Patch did not divide into halves based on topology and an "
157  << "included angle of " << includedAngle << " degrees"
158  << exit(FatalError);
159  }
160 
161  // Convert the zone labels into bools describing which half each face is in
162  boolList faceHalves(pp.size());
163  forAll(faceHalves, facei)
164  {
165  faceHalves[facei] = faceZones[facei] == 1;
166  }
167 
168  return faceHalves;
169 }
170 
171 } // End namespace Foam
172 
173 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
174 
175 void Foam::polyMeshUnMergeCyclics(polyMesh& mesh, const scalar includedAngle)
176 {
177  // Add cyclic patches. Clone all the existing patches into a new list, and
178  // add two additional empty cyclic patches for each merged-cyclic. Then
179  // re-patch the mesh
181 
182  forAll(mesh.boundaryMesh(), patchi)
183  {
184  const polyPatch& pp = mesh.boundaryMesh()[patchi];
185  patches.append
186  (
187  pp.clone
188  (
189  mesh.boundaryMesh(),
190  patches.size(),
191  pp.size(),
192  pp.start()
193  ).ptr()
194  );
195  }
196 
197  bool hasMergedCyclics = false;
198  labelList patchHalf0(mesh.boundaryMesh().size(), -1);
199  labelList patchHalf1(mesh.boundaryMesh().size(), -1);
200 
201  forAll(mesh.boundaryMesh(), patchi)
202  {
203  const polyPatch& pp = mesh.boundaryMesh()[patchi];
204  if (isA<mergedCyclicPolyPatch>(pp))
205  {
206  hasMergedCyclics = true;
207 
208  patchHalf0[patchi] = patches.size();
209 
210  patches.append
211  (
212  new cyclicPolyPatch
213  (
214  pp.name() + "-half-0",
215  0,
216  mesh.nFaces(),
217  patches.size(),
218  mesh.boundaryMesh(),
219  cyclicPolyPatch::typeName,
220  pp.name() + "-half-1"
221  )
222  );
223 
224  patchHalf1[patchi] = patches.size();
225 
226  patches.append
227  (
228  new cyclicPolyPatch
229  (
230  pp.name() + "-half-1",
231  0,
232  mesh.nFaces(),
233  patches.size(),
234  mesh.boundaryMesh(),
235  cyclicPolyPatch::typeName,
236  pp.name() + "-half-0"
237  )
238  );
239  }
240  }
241 
242  if (!hasMergedCyclics)
243  {
244  return;
245  }
246 
247  mesh.removeBoundary();
248  mesh.addPatches(patches);
249 
250  // Move faces from the merged cyclics to the pairs of cyclic patches
251  polyTopoChange meshMod(mesh, false);
252  forAll(mesh.boundaryMesh(), patchi)
253  {
254  const polyPatch& pp = mesh.boundaryMesh()[patchi];
255  if (isA<mergedCyclicPolyPatch>(pp))
256  {
257  const boolList patchFaceHalves
258  (
259  primitivePatchGetHalves(pp, includedAngle)
260  );
261 
262  forAll(pp, patchFacei)
263  {
264  const label facei = pp.start() + patchFacei;
265 
266  meshMod.modifyFace
267  (
268  mesh.faces()[facei],
269  facei,
270  mesh.faceOwner()[facei],
271  -1,
272  false,
273  patchFaceHalves[patchFacei]
274  ? patchHalf0[patchi]
275  : patchHalf1[patchi],
276  -1,
277  false
278  );
279  }
280  }
281  }
282  autoPtr<polyTopoChangeMap> map = meshMod.changeMesh(mesh, false);
283 
284  // Remove the (now) empty merged cyclic patches. Copy everything except the
285  // merged cyclics into a patch list and use this then re-patch the mesh.
286  patches.clear();
287 
288  forAll(mesh.boundaryMesh(), patchi)
289  {
290  const polyPatch& pp = mesh.boundaryMesh()[patchi];
291  if (!isA<mergedCyclicPolyPatch>(pp))
292  {
293  patches.append
294  (
295  pp.clone
296  (
297  mesh.boundaryMesh(),
298  patches.size(),
299  pp.size(),
300  pp.start()
301  ).ptr()
302  );
303  }
304  }
305 
306  mesh.removeBoundary();
307  mesh.addPatches(patches);
308 }
309 
310 
311 // ************************************************************************* //
const fvPatchList & patches
const polyBoundaryMesh & boundaryMesh() const
Return boundary mesh.
Definition: polyMesh.H:453
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
autoPtr< polyTopoChangeMap > changeMesh(polyMesh &mesh, const bool inflate, const bool syncParallel=true, const bool orderCells=false, const bool orderPoints=false)
Inplace changes mesh without change of patches.
const word & name() const
Return name.
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:124
error FatalError
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:306
labelList primitivePatchGetZones(const primitivePatch &pp, const scalar includedAngle)
label nFaces() const
Unit conversion functions.
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:164
static int compare(const edge &, const edge &)
Compare edges.
Definition: edgeI.H:36
label fcIndex(const label i) const
Return the forward circular index, i.e. the next index.
Definition: UListI.H:58
label findMin(const ListType &, const label start=0)
Find index of min element (and less than given element).
fvMesh & mesh
scalar degToRad(const scalar deg)
Conversion from degrees to radians.
void polyMeshUnMergeCyclics(polyMesh &mesh, const scalar includedAngle=165)
Find all patches of type "mergedCyclic" in the given mesh and split them.
virtual autoPtr< polyPatch > clone(const polyBoundaryMesh &bm) const
Construct and return a clone, resetting the boundary mesh.
Definition: polyPatch.H:225
A list of faces which address into the list of points.
dimensionedScalar cos(const dimensionedScalar &ds)
A 1D vector of objects of type <T> that resizes itself as necessary to accept the new objects...
Definition: DynamicList.H:56
Cyclic plane patch.
static int compare(const face &, const face &)
Compare faces.
Definition: face.C:48
virtual const labelList & faceOwner() const
Return face owner.
Definition: polyMesh.C:1237
DynamicList< T, SizeInc, SizeMult, SizeDiv > & append(const T &)
Append an element at the end of the list.
Definition: DynamicListI.H:296
const labelListList & edgeFaces() const
Return edge-face addressing.
List< label > labelList
A List of labels.
Definition: labelList.H:56
virtual const faceList & faces() const
Return raw faces.
Definition: polyMesh.C:1224
const Field< PointType > & faceNormals() const
Return face normals for patch.
boolList primitivePatchGetHalves(const primitivePatch &pp, const scalar includedAngle)
void addPatches(const List< polyPatch *> &, const bool validBoundary=true)
Add boundary patches.
Definition: polyMesh.C:989
label findMax(const ListType &, const label start=0)
Find index of max element (and larger than given element).
void removeBoundary()
Remove boundary patches.
Definition: polyMeshClear.C:36
label findIndex(const ListType &, typename ListType::const_reference, const label start=0)
Find first occurrence of given element and return index,.
label size() const
Return the number of elements in the UPtrList.
Definition: UPtrListI.H:29
void modifyFace(const face &f, const label facei, const label own, const label nei, const bool flipFaceFlux, const label patchID, const label zoneID, const bool zoneFlip)
Modify vertices or cell of face.
label patchi
Direct mesh changes based on v1.3 polyTopoChange syntax.
label start() const
Return start label of this patch in the polyMesh face list.
Definition: polyPatch.H:306
const labelListList & faceEdges() const
Return face-edge addressing.
An auto-pointer similar to the STL auto_ptr but with automatic casting to a reference to the type and...
Definition: PtrList.H:52
Mesh consisting of general polyhedral cells.
Definition: polyMesh.H:76
A patch is a list of labels that address the faces in the global face list.
Definition: polyPatch.H:66
void clear()
Clear the addressed list, i.e. set the size to zero.
Definition: DynamicListI.H:236
Namespace for OpenFOAM.