cutPoly.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) 2022-2023 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 "cutPolyValue.H"
27 #include "OBJstream.H"
28 
29 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
30 
31 namespace Foam
32 {
33 
34 template<class Type>
35 Type min(const UIndirectList<Type>& l)
36 {
37  Type result = pTraits<Type>::max;
38  forAll(l, i)
39  {
40  result = min(result, l[i]);
41  }
42  return result;
43 }
44 
45 template<class Type>
46 Type max(const UIndirectList<Type>& l)
47 {
48  Type result = pTraits<Type>::min;
49  forAll(l, i)
50  {
51  result = max(result, l[i]);
52  }
53  return result;
54 }
55 
56 }
57 
58 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
59 
61 (
62  const face& f,
63  const scalarField& pAlphas,
64  const scalar isoAlpha
65 )
66 {
67  UIndirectList<scalar> fAlphas(pAlphas, f);
68 
69  // Quick reject if all alpha values are above or below the iso value
70  if (min(fAlphas) > isoAlpha || isoAlpha > max(fAlphas))
71  {
72  return List<labelPair>();
73  }
74 
75  // Find the starting point
76  label fpi0 = 0;
77  while ((fAlphas[fpi0] < isoAlpha) == separatedBelow)
78  {
79  ++ fpi0;
80  }
81 
82  // Create cuts on every edge for which there is a sign change
83  DynamicList<labelPair> cuts(1);
84  label cuti = 0;
85  forAll(f, i)
86  {
87  const label fpi1 = f.fcIndex(fpi0);
88 
89  if ((fAlphas[fpi0] < isoAlpha) != (fAlphas[fpi1] < isoAlpha))
90  {
91  if (cuti == 0)
92  {
93  cuts.append({fpi0, -1});
94  }
95  else
96  {
97  cuts.last()[1] = fpi0;
98  }
99 
100  cuti = cuts.last().fcIndex(cuti);
101  }
102 
103  fpi0 = fpi1;
104  }
105 
106  // There should have been an even number of cuts
107  if (cuti != 0)
108  {
110  << "Cutting values " << fAlphas << " with iso-value " << isoAlpha
111  << " resulted in an odd number of cuts"
112  << exit(FatalError);
113  }
114 
115  cuts.shrink();
116  return cuts;
117 }
118 
119 
121 (
122  const cell& c,
123  const cellEdgeAddressing& cAddr,
124  const faceUList& fs,
125  const List<List<labelPair>>& fCuts,
126  const scalarField& pAlphas,
127  const scalar isoAlpha
128 )
129 {
130  // Quick return if not cut
131  label nCellFaceCuts = 0;
132  forAll(c, cfi)
133  {
134  nCellFaceCuts += fCuts[c[cfi]].size();
135  }
136  if (nCellFaceCuts == 0)
137  {
138  return labelListList();
139  }
140 
141  // Get local cell-face-edge addressing
142  const CompactListList<label>& cfiAndFeiToCei = cAddr.cfiAndFeiToCei();
143  const List<Pair<labelPair>>& ceiToCfiAndFei = cAddr.ceiToCfiAndFei();
144  const boolList& cOwns = cAddr.cOwns();
145 
146  // For each cell-face and face-edge, get the face-edge at the other end of
147  // the edge's cut, or -1 if the edge is not cut. This allows us to walk
148  // along face cuts. It also doubles as the "visited" array during the walk.
149  // As cuts are visited, the relevant label in this array is set to -1 to
150  // ensure that the cut is not considered twice.
151  labelListList cellFaceAndFaceEdgeCutToFaceEdge(c.size());
152  forAll(c, cfi)
153  {
154  cellFaceAndFaceEdgeCutToFaceEdge[cfi] =
155  labelList(fs[c[cfi]].size(), -1);
156 
157  forAll(fCuts[c[cfi]], fci)
158  {
159  const label fei0 = fCuts[c[cfi]][fci].first();
160  const label fei1 = fCuts[c[cfi]][fci].second();
161 
162  if (cOwns[cfi])
163  {
164  cellFaceAndFaceEdgeCutToFaceEdge[cfi][fei0] = fei1;
165  }
166  else
167  {
168  cellFaceAndFaceEdgeCutToFaceEdge[cfi][fei1] = fei0;
169  }
170  }
171  }
172 
173  // Walk around the face cuts to generate the cell cuts
175  {
176  label nCellFaceCutsVisited = 0;
177  label cfi0 = 0, fei0 = 0;
178  while (nCellFaceCutsVisited < nCellFaceCuts)
179  {
180  // Find the next unvisited connection
181  bool found = false;
182  for (label cfj0 = cfi0; cfj0 < c.size(); ++ cfj0)
183  {
184  for (label fej0 = fei0; fej0 < fs[c[cfj0]].size(); ++ fej0)
185  {
186  const label fej1 =
187  cellFaceAndFaceEdgeCutToFaceEdge[cfj0][fej0];
188 
189  if (fej1 != -1)
190  {
191  cfi0 = cfj0;
192  fei0 = fej0;
193  found = true;
194  }
195 
196  if (found) break;
197  }
198 
199  if (found) break;
200 
201  fei0 = 0;
202  }
203  if (!found)
204  {
206  << "Could not find next unvisited connection for cell cut"
207  << exit(FatalError);
208  }
209 
210  // Walk around the cell from face to face to form a cut
211  label cfi = cfi0, fei = fei0;
212  DynamicList<label> cut(8);
213  while (cellFaceAndFaceEdgeCutToFaceEdge[cfi][fei] != -1)
214  {
215  ++ nCellFaceCutsVisited;
216 
217  const label otherFei =
218  cellFaceAndFaceEdgeCutToFaceEdge[cfi][fei];
219 
220  cellFaceAndFaceEdgeCutToFaceEdge[cfi][fei] = -1;
221 
222  const label cei = cfiAndFeiToCei[cfi][otherFei];
223 
224  cut.append(cei);
225 
226  const labelPair next =
227  ceiToCfiAndFei[cei].other({cfi, otherFei});
228 
229  cfi = next.first();
230  fei = next.second();
231  }
232 
233  // Copy the data over into the result list
234  cuts.append(labelList());
235  cuts.last().transfer(cut);
236  }
237  }
238 
239  cuts.shrink();
240  return cuts;
241 }
242 
243 
245 (
246  const face& f,
247  const List<labelPair>& fCuts,
248  const pointField& ps,
249  const scalarField& pAlphas,
250  const scalar isoAlpha,
251  OBJstream& obj
252 )
253 {
254  forAll(fCuts, cuti)
255  {
256  pointField cutPs(2);
257 
258  forAll(fCuts[cuti], i)
259  {
260  const edge e = f.faceEdge(fCuts[cuti][i]);
261 
262  cutPs[i] = edgeCutValue(e, pAlphas, isoAlpha, ps);
263  }
264 
265  obj.write(edge(0, 1), cutPs);
266  }
267 }
268 
269 
271 (
272  const cell& c,
273  const cellEdgeAddressing& cAddr,
274  const List<List<label>>& cCuts,
275  const faceList& fs,
276  const pointField& ps,
277  const scalarField& pAlphas,
278  const scalar isoAlpha,
279  OBJstream& obj
280 )
281 {
282  // Quick return if not cut
283  if (cCuts.size() == 0)
284  {
285  return;
286  }
287 
288  // Get local cell-face-edge addressing
289  const List<Pair<labelPair>>& ceiToCfiAndFei = cAddr.ceiToCfiAndFei();
290 
291  // Write out each cut as a single face
292  forAll(cCuts, cuti)
293  {
294  if (cCuts[cuti].size() < 3) continue;
295 
296  pointField cutPs(cCuts[cuti].size());
297 
298  forAll(cCuts[cuti], i)
299  {
300  const label cei = cCuts[cuti][i];
301  const label cfi = ceiToCfiAndFei[cei][0][0];
302  const label fei = ceiToCfiAndFei[cei][0][1];
303 
304  const edge e = fs[c[cfi]].faceEdge(fei);
305 
306  cutPs[i] = edgeCutValue(e, pAlphas, isoAlpha, ps);
307  }
308 
309  obj.write(face(identityMap(cCuts[cuti].size())), cutPs, false);
310  }
311 }
312 
313 
314 // ************************************************************************* //
bool found
#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
DynamicList< T, SizeInc, SizeMult, SizeDiv > & append(const T &)
Append an element at the end of the list.
Definition: DynamicListI.H:296
DynamicList< T, SizeInc, SizeMult, SizeDiv > & shrink()
Shrink the allocated space to the number of elements used.
Definition: DynamicListI.H:252
A 1D array of objects of type <T>, where the size of the vector is known and used for subscript bound...
Definition: List.H:91
OFstream which keeps track of vertices.
Definition: OBJstream.H:56
virtual Ostream & write(const char)
Write character.
Definition: OBJstream.C:81
const Type & second() const
Return second.
Definition: Pair.H:110
const Type & first() const
Return first.
Definition: Pair.H:98
A List with indirect addressing.
Definition: UIndirectList.H:60
A 1D vector of objects of type <T>, where the size of the vector is known and can be used for subscri...
Definition: UList.H:74
label size() const
Return the number of elements in the UList.
Definition: UListI.H:311
T & last()
Return the last element of the list.
Definition: UListI.H:128
const List< Pair< labelPair > > & ceiToCfiAndFei() const
Map from cell-edge-index to cell-face-index and face-edge-index.
const boolList & cOwns() const
For each cell-face, whether or not the cell owns it.
const CompactListList< label > & cfiAndFeiToCei() const
Map from cell-face-index and face-edge-index to cell-edge-index.
A cell is defined as a list of faces with extra functionality.
Definition: cell.H:60
An edge is a list of two point labels. The functionality it provides supports the discretisation on a...
Definition: edge.H:61
A face is a list of labels corresponding to mesh vertices.
Definition: face.H:76
Traits class for primitives.
Definition: pTraits.H:53
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:334
const dimensionedScalar e
Elementary charge.
const dimensionedScalar c
Speed of light in a vacuum.
Type edgeCutValue(const edge &e, const scalar lambda, const Field< Type > &pPsis)
Linearly interpolate a value from the end points to the cut point of an.
static const bool separatedBelow
This flag determines which side of the iso-surface is separated into.
Definition: cutPoly.H:65
void writeCellCuts(const cell &c, const cellEdgeAddressing &cAddr, const List< List< label >> &cCuts, const faceList &fs, const pointField &ps, const scalarField &pAlphas, const scalar isoAlpha, OBJstream &obj)
Write the cuts for a given cell to an OBJ file.
Definition: cutPoly.C:271
void writeFaceCuts(const face &f, const List< labelPair > &fCuts, const pointField &ps, const scalarField &pAlphas, const scalar isoAlpha, OBJstream &obj)
Write the cuts for a given face to an OBJ file.
Definition: cutPoly.C:245
List< labelPair > faceCuts(const face &f, const scalarField &pAlphas, const scalar isoAlpha)
Return the cuts for a given face. This returns a list of pairs of.
Definition: cutPoly.C:61
labelListList cellCuts(const cell &c, const cellEdgeAddressing &cAddr, const faceUList &fs, const List< List< labelPair >> &fCuts, const scalarField &pAlphas, const scalar isoAlpha)
Return the cuts for a given cell. This returns a list of lists of cell-edge.
Definition: cutPoly.C:121
Namespace for OpenFOAM.
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:124
List< label > labelList
A List of labels.
Definition: labelList.H:56
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
Type max(const UIndirectList< Type > &l)
Definition: cutPoly.C:46
layerAndWeight min(const layerAndWeight &a, const layerAndWeight &b)
List< labelList > labelListList
A List of labelList.
Definition: labelList.H:57
layerAndWeight max(const layerAndWeight &a, const layerAndWeight &b)
error FatalError
labelList identityMap(const label len)
Create identity map (map[i] == i) of given length.
Definition: ListOps.C:104
Type min(const UIndirectList< Type > &l)
Definition: cutPoly.C:35
labelList f(nPoints)