isoSurfaceBase.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 | www.openfoam.com
6 \\/ M anipulation |
7-------------------------------------------------------------------------------
8 Copyright (C) 2019-2020 OpenCFD Ltd.
9-------------------------------------------------------------------------------
10License
11 This file is part of OpenFOAM.
12
13 OpenFOAM is free software: you can redistribute it and/or modify it
14 under the terms of the GNU General Public License as published by
15 the Free Software Foundation, either version 3 of the License, or
16 (at your option) any later version.
17
18 OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
19 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
20 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
21 for more details.
22
23 You should have received a copy of the GNU General Public License
24 along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
25
26\*---------------------------------------------------------------------------*/
27
28#include "isoSurfaceBase.H"
29#include "polyMesh.H"
30#include "tetMatcher.H"
31#include "cyclicACMIPolyPatch.H"
32
33// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
34
37
38
39// * * * * * * * * * * * * * * * Local Functions * * * * * * * * * * * * * * //
40
41namespace Foam
42{
43
44// Test face for edge cuts
45inline static bool isFaceCut
46(
47 const scalar isoval,
48 const scalarField& pointValues,
49 const labelUList& indices
50)
51{
52 auto iter = indices.cbegin();
53 const auto last = indices.cend();
54
55 // if (iter == last) return false; // ie, indices.empty()
56
57 const bool aLower = (pointValues[*iter] < isoval);
58 ++iter;
59
60 while (iter != last)
61 {
62 if (aLower != (pointValues[*iter] < isoval))
63 {
64 return true;
65 }
66 ++iter;
67 }
68
69 return false;
70}
71
72} // End namespace Foam
73
74
75// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
76
78(
79 const polyMesh& mesh,
80 const scalarField& cellValues,
81 const scalarField& pointValues,
82 const scalar iso,
83 const isoSurfaceParams& params
84)
85:
87 isoSurfaceParams(params),
88 mesh_(mesh),
89 cVals_(cellValues),
90 pVals_(pointValues),
91 iso_(iso),
92 ignoreBoundaryFaces_(),
93 meshCells_()
94{}
95
96
97// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
98
100{
101 // Determine boundary pyramids to ignore (originating from ACMI faces)
102 // Maybe needs to be optional argument or more general detect colocated
103 // faces.
104
105 for (const polyPatch& pp : mesh_.boundaryMesh())
106 {
107 if (isA<cyclicACMIPolyPatch>(pp))
108 {
109 ignoreBoundaryFaces_.set(labelRange(pp.offset(), pp.size()));
110 }
111 }
112}
113
114
116(
117 const UList<cutType>& cuts,
118 const uint8_t maskValue
119)
120{
121 label count = 0;
122
123 for (const cutType cut : cuts)
124 {
125 if (maskValue ? (cut & maskValue) != 0 : !cut)
126 {
127 ++count;
128 }
129 }
130
131 return count;
132}
133
134
136{
137 for (cutType& cut : cuts)
138 {
139 if (cut != cutType::BLOCKED)
140 {
141 cut = cutType::UNVISITED;
142 }
143 }
144}
145
146
148(
149 UList<cutType>& cuts,
150 const bitSet& ignoreCells
151) const
152{
153 label count = 0;
154
155 for (const label celli : ignoreCells)
156 {
157 if (celli >= cuts.size())
158 {
159 break;
160 }
161
162 cuts[celli] = cutType::BLOCKED;
163 ++count;
164 }
165
166 return count;
167}
168
169
171(
172 UList<cutType>& cuts,
173 const boundBox& bb,
174 const volumeType::type volType
175) const
176{
177 label count = 0;
178
179 // Mark cells inside/outside bounding box as blocked
180 const bool keepInside = (volType == volumeType::INSIDE);
181
182 if (!keepInside && volType != volumeType::OUTSIDE)
183 {
184 // Could warn about invalid...
185 }
186 else if (bb.valid())
187 {
188 const pointField& cc = mesh_.cellCentres();
189
190 forAll(cuts, celli)
191 {
192 if
193 (
194 cuts[celli] == cutType::UNVISITED
195 && (bb.contains(cc[celli]) ? keepInside : !keepInside)
196 )
197 {
198 cuts[celli] = cutType::BLOCKED;
199 ++count;
200 }
201 }
202 }
203
204 return count;
205}
206
207
209{
210 // Don't consider SPHERE cuts in the total number of cells cut
211 constexpr uint8_t realCut(cutType::CUT | cutType::TETCUT);
212
213 cuts.resize(mesh_.nCells(), cutType::UNVISITED);
214
215 label nCuts = 0;
216 forAll(cuts, celli)
217 {
218 if (cuts[celli] == cutType::UNVISITED)
219 {
220 cuts[celli] = getCellCutType(celli);
221
222 if ((cuts[celli] & realCut) != 0)
223 {
224 ++nCuts;
225 }
226 }
227 }
228
229 return nCuts;
230}
231
232
235{
236 return
237 (
238 (
239 mesh_.isInternalFace(facei)
240 || !ignoreBoundaryFaces_.test(facei-mesh_.nInternalFaces())
241 )
242 && isFaceCut(iso_, pVals_, mesh_.faces()[facei])
243 ) ? cutType::CUT : cutType::NOTCUT;
244}
245
246
249{
250 // Tet version
251 if (tetMatcher::test(mesh_, celli))
252 {
253 for (const label facei : mesh_.cells()[celli])
254 {
255 if
256 (
257 !mesh_.isInternalFace(facei)
258 && ignoreBoundaryFaces_.test(facei-mesh_.nInternalFaces())
259 )
260 {
261 continue;
262 }
263
264 if (isFaceCut(iso_, pVals_, mesh_.faces()[facei]))
265 {
266 return cutType::TETCUT;
267 }
268 }
269
270 return cutType::NOTCUT;
271 }
272
273
274 // Regular cell
275 label nPyrEdges = 0;
276 label nPyrCuts = 0;
277
278 const bool cellLower = (cVals_[celli] < iso_);
279
280 for (const label facei : mesh_.cells()[celli])
281 {
282 if
283 (
284 !mesh_.isInternalFace(facei)
285 && ignoreBoundaryFaces_.test(facei-mesh_.nInternalFaces())
286 )
287 {
288 continue;
289 }
290
291 const face& f = mesh_.faces()[facei];
292
293 // Count pyramid edges cut
294 for (const label pointi : f)
295 {
296 ++nPyrEdges;
297
298 if (cellLower != (pVals_[pointi] < iso_))
299 {
300 ++nPyrCuts;
301 }
302 }
303 }
304
305 if (nPyrCuts == 0)
306 {
307 return cutType::NOTCUT;
308 }
309
310 // If all pyramid edges are cut (no outside faces),
311 // it is a sphere cut
312
313 return (nPyrCuts == nPyrEdges) ? cutType::SPHERE : cutType::CUT;
314}
315
316
317// ************************************************************************* //
A 1D array of objects of type <T>, where the size of the vector is known and used for subscript bound...
Definition: List.H:77
void resize(const label len)
Adjust allocated size of list.
Definition: ListI.H:139
const_iterator cend() const noexcept
Return const_iterator to end traversing the constant UList.
Definition: UListI.H:364
const_iterator cbegin() const noexcept
Return const_iterator to begin traversing the constant UList.
Definition: UListI.H:343
void size(const label n)
Older name for setAddressableSize.
Definition: UList.H:114
A bitSet stores bits (elements with only two states) in packed internal format and supports a variety...
Definition: bitSet.H:66
A bounding box defined in terms of min/max extrema points.
Definition: boundBox.H:64
bool valid() const
Bounding box is non-inverted.
Definition: boundBoxI.H:76
bool contains(const point &pt) const
Contains point? (inside or on edge)
Definition: boundBoxI.H:271
A face is a list of labels corresponding to mesh vertices.
Definition: face.H:75
labelListList blockCells() const
Return block cells.
Low-level components common to various iso-surface algorithms.
static void resetCuts(UList< cutType > &cuts)
Restore non-BLOCKED state to an UNVISITED state.
static label countCutType(const UList< cutType > &cuts, const uint8_t maskValue)
Count the number of cuts matching the mask type.
cutType
The type of cell/face cuts.
label calcCellCuts(List< cutType > &cuts) const
Populate a list of candidate cell cuts using getCellCutType()
cutType getFaceCutType(const label facei) const
Determine face cut for an individual face.
void ignoreCyclics()
Set ignoreBoundaryFaces to ignore cyclics (cyclicACMI)
cutType getCellCutType(const label celli) const
Preferences for controlling iso-surface algorithms.
A range or interval of labels defined by a start and a size.
Definition: labelRange.H:58
Mesh consisting of general polyhedral cells.
Definition: polyMesh.H:81
A patch is a list of labels that address the faces in the global face list.
Definition: polyPatch.H:75
const polyBoundaryMesh & boundaryMesh() const
Return boundaryMesh reference.
Definition: polyPatch.C:315
Patchify triangles based on orientation w.r.t other (triangulated or triangulatable) surfaces.
static bool test(const UList< face > &faces)
Definition: tetMatcher.C:89
type
Volume classification types.
Definition: volumeType.H:66
@ OUTSIDE
A location outside the volume.
Definition: volumeType.H:69
@ INSIDE
A location inside the volume.
Definition: volumeType.H:68
dynamicFvMesh & mesh
Convenience macros for instantiating iso-surface interpolate methods.
#define defineIsoSurfaceInterpolateMethods(ThisClass)
Namespace for OpenFOAM.
static bool isFaceCut(const scalar isoval, const scalarField &pointValues, const labelUList &indices)
labelList f(nPoints)
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:333