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-2024 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 " << radToDeg(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  );
277  }
278  }
279  }
280  autoPtr<polyTopoChangeMap> map = meshMod.changeMesh(mesh);
281 
282  // Remove the (now) empty merged cyclic patches. Copy everything except the
283  // merged cyclics into a patch list and use this then re-patch the mesh.
284  patches.clear();
285 
286  forAll(mesh.boundaryMesh(), patchi)
287  {
288  const polyPatch& pp = mesh.boundaryMesh()[patchi];
289  if (!isA<mergedCyclicPolyPatch>(pp))
290  {
291  patches.append
292  (
293  pp.clone
294  (
295  mesh.boundaryMesh(),
296  patches.size(),
297  pp.size(),
298  pp.start()
299  ).ptr()
300  );
301  }
302  }
303 
304  mesh.removeBoundary();
305  mesh.addPatches(patches);
306 }
307 
308 
309 // ************************************************************************* //
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
A 1D vector of objects of type <T> that resizes itself as necessary to accept the new objects.
Definition: DynamicList.H:78
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:164
const labelListList & edgeFaces() const
Return edge-face addressing.
const labelListList & faceEdges() const
Return face-edge addressing.
const Field< PointType > & faceNormals() const
Return face normals for patch.
label size() const
Return the number of elements in the UList.
Definition: UListI.H:311
label fcIndex(const label i) const
Return the forward circular index, i.e. the next index.
Definition: UListI.H:58
label size() const
Return the number of elements in the UPtrList.
Definition: UPtrListI.H:29
An auto-pointer similar to the STL auto_ptr but with automatic casting to a reference to the type and...
Definition: autoPtr.H:51
Cyclic plane patch.
static int compare(const edge &, const edge &)
Compare edges.
Definition: edgeI.H:36
static int compare(const face &, const face &)
Compare faces.
Definition: face.C:48
const word & name() const
Return name.
Mesh consisting of general polyhedral cells.
Definition: polyMesh.H:80
virtual const faceList & faces() const
Return raw faces.
Definition: polyMesh.C:1326
virtual const labelList & faceOwner() const
Return face owner.
Definition: polyMesh.C:1339
const polyBoundaryMesh & boundaryMesh() const
Return boundary mesh.
Definition: polyMesh.H:404
void addPatches(const List< polyPatch * > &, const bool validBoundary=true)
Add boundary patches.
Definition: polyMesh.C:1084
void removeBoundary()
Remove boundary patches.
Definition: polyMeshClear.C:36
A patch is a list of labels that address the faces in the global face list.
Definition: polyPatch.H:70
label start() const
Return start label of this patch in the polyMesh face list.
Definition: polyPatch.H:280
virtual autoPtr< polyPatch > clone(const polyBoundaryMesh &bm) const
Construct and return a clone, resetting the boundary mesh.
Definition: polyPatch.H:215
Direct mesh changes based on v1.3 polyTopoChange syntax.
autoPtr< polyTopoChangeMap > changeMesh(polyMesh &mesh, const bool syncParallel=true, const bool orderCells=false, const bool orderPoints=false)
Inplace changes mesh without change of patches.
void modifyFace(const face &f, const label facei, const label own, const label nei, const bool flipFaceFlux, const label patchID)
Modify vertices or cell of face.
label nFaces() const
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:334
label patchi
const fvPatchList & patches
Namespace for OpenFOAM.
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:124
scalar radToDeg(const scalar rad)
Convert radians to degrees.
scalar degToRad(const scalar deg)
Convert degrees to radians.
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
void polyMeshUnMergeCyclics(polyMesh &mesh, const scalar includedAngle=degToRad(165))
Find all patches of type "mergedCyclic" in the given mesh and split them.
label findMax(const ListType &, const label start=0)
Find index of max element (and larger than given element).
boolList primitivePatchGetHalves(const primitivePatch &pp, const scalar includedAngle)
labelList primitivePatchGetZones(const primitivePatch &pp, const scalar includedAngle)
label findMin(const ListType &, const label start=0)
Find index of min element (and less than given element).
label findIndex(const ListType &, typename ListType::const_reference, const label start=0)
Find first occurrence of given element and return index,.
error FatalError
dimensionedScalar cos(const dimensionedScalar &ds)