isoSurface.H 11.5 KB
Newer Older
mattijs's avatar
mattijs committed
1
2
3
4
/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
Mark Olesen's avatar
Mark Olesen committed
5
    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
mattijs's avatar
mattijs committed
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
     \\/     M anipulation  |
-------------------------------------------------------------------------------
License
    This file is part of OpenFOAM.

    OpenFOAM is free software; you can redistribute it and/or modify it
    under the terms of the GNU General Public License as published by the
    Free Software Foundation; either version 2 of the License, or (at your
    option) any later version.

    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
    for more details.

    You should have received a copy of the GNU General Public License
    along with OpenFOAM; if not, write to the Free Software Foundation,
    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA

Class
    Foam::isoSurface

Description
    A surface formed by the iso value.
mattijs's avatar
mattijs committed
30
    After "Regularised Marching Tetrahedra: improved iso-surface extraction",
mattijs's avatar
mattijs committed
31
32
33
34
35
36
37
38
39
    G.M. Treece, R.W. Prager and A.H. Gee.

    Note:
    - in parallel the regularisation (coarsening) always takes place
    and slightly different surfaces will be created compared to non-parallel.
    The surface will still be continuous though!
    - does tets without using cell centres/cell values. Not tested.
    - regularisation can create duplicate triangles/non-manifold surfaces.
    Current handling of those is bit ad-hoc for now and not perfect.
mattijs's avatar
mattijs committed
40
41
42
    - uses geometric merge with fraction of bounding box as distance.
    - triangles can be between two cell centres so constant sampling
      does not make sense.
mattijs's avatar
mattijs committed
43
    - on empty patches behaves like zero gradient.
mattijs's avatar
mattijs committed
44
    - does not do 2D correctly, creates non-flat iso surface.
mattijs's avatar
mattijs committed
45
46
47
48
49
50
51
52
53
54

SourceFiles
    isoSurface.C

\*---------------------------------------------------------------------------*/

#ifndef isoSurface_H
#define isoSurface_H

#include "triSurface.H"
mattijs's avatar
mattijs committed
55
56
#include "labelPair.H"
#include "pointIndexHit.H"
57
#include "PackedBoolList.H"
mattijs's avatar
mattijs committed
58
#include "volFields.H"
mattijs's avatar
mattijs committed
59
60
61
62
63
64

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

namespace Foam
{

mattijs's avatar
mattijs committed
65
class fvMesh;
mattijs's avatar
mattijs committed
66
67
68
69
70
71
72
73
74
75
76

/*---------------------------------------------------------------------------*\
                       Class isoSurface Declaration
\*---------------------------------------------------------------------------*/

class isoSurface
:
    public triSurface
{
    // Private data

mattijs's avatar
mattijs committed
77
78
79
80
81
82
83
84
85
86
87
88
89
        enum segmentCutType
        {
            NEARFIRST,      // intersection close to e.first()
            NEARSECOND,     //  ,,                   e.second()
            NOTNEAR         // no intersection
        };

        enum cellCutType
        {
            NOTCUT,     // not cut
            SPHERE,     // all edges to cell centre cut
            CUT         // normal cut
        };
mattijs's avatar
mattijs committed
90
91


mattijs's avatar
mattijs committed
92
        //- Reference to mesh
mattijs's avatar
mattijs committed
93
        const fvMesh& mesh_;
mattijs's avatar
mattijs committed
94
95
96
97

        //- Isosurface value
        const scalar iso_;

mattijs's avatar
mattijs committed
98
99
100
        //- When to merge points
        const scalar mergeDistance_;

mattijs's avatar
mattijs committed
101

mattijs's avatar
mattijs committed
102
103
104
105
106
107
108
109
110
        //- Whether face might be cut
        List<cellCutType> faceCutType_;

        //- Whether cell might be cut
        List<cellCutType> cellCutType_;

        //- Estimated number of cut cells
        label nCutCells_;

mattijs's avatar
mattijs committed
111
112
113
114
115
116
117
118
119
        //- For every triangle the original cell in mesh
        labelList meshCells_;

        //- For every unmerged triangle point the point in the triSurface
        labelList triPointMergeMap_;


    // Private Member Functions

mattijs's avatar
mattijs committed
120
121
        //- Get location of iso value as fraction inbetween s0,s1
        scalar isoFraction
mattijs's avatar
mattijs committed
122
123
124
        (
            const scalar s0,
            const scalar s1
mattijs's avatar
mattijs committed
125
126
        ) const;

mattijs's avatar
mattijs committed
127
128
129
130
131
132
133
134
135
        //- Check if any edge of a face is cut
        bool isEdgeOfFaceCut
        (
            const scalarField& pVals,
            const face& f,
            const bool ownLower,
            const bool neiLower
        ) const;

mattijs's avatar
mattijs committed
136
137
138
139
140
141
142
143
144
145
        void getNeighbour
        (
            const labelList& boundaryRegion,
            const volScalarField& cVals,
            const label cellI,
            const label faceI,
            scalar& nbrValue,
            point& nbrPoint
        ) const;

mattijs's avatar
mattijs committed
146
147
        //- Set faceCutType,cellCutType.
        void calcCutTypes
mattijs's avatar
mattijs committed
148
        (
mattijs's avatar
mattijs committed
149
            const labelList& boundaryRegion,
mattijs's avatar
mattijs committed
150
151
152
            const volScalarField& cVals,
            const scalarField& pVals
        );
mattijs's avatar
mattijs committed
153
154
155
156
157
158
159
160
161
162
163
164
165

        static labelPair findCommonPoints
        (
            const labelledTri&,
            const labelledTri&
        );

        static point calcCentre(const triSurface&);

        static pointIndexHit collapseSurface
        (
            pointField& localPoints,
            DynamicList<labelledTri, 64>& localTris
mattijs's avatar
mattijs committed
166
167
        );

mattijs's avatar
mattijs committed
168
169
170
        //- Determine per cc whether all near cuts can be snapped to single
        //  point.
        void calcSnappedCc
mattijs's avatar
mattijs committed
171
        (
mattijs's avatar
mattijs committed
172
173
            const labelList& boundaryRegion,
            const volScalarField& cVals,
mattijs's avatar
mattijs committed
174
            const scalarField& pVals,
mattijs's avatar
mattijs committed
175
            DynamicList<point>& snappedPoints,
mattijs's avatar
mattijs committed
176
177
178
179
180
181
182
            labelList& snappedCc
        ) const;

        //- Determine per point whether all near cuts can be snapped to single
        //  point.
        void calcSnappedPoint
        (
183
            const PackedBoolList& isBoundaryPoint,
mattijs's avatar
mattijs committed
184
185
            const labelList& boundaryRegion,
            const volScalarField& cVals,
mattijs's avatar
mattijs committed
186
            const scalarField& pVals,
mattijs's avatar
mattijs committed
187
            DynamicList<point>& snappedPoints,
mattijs's avatar
mattijs committed
188
189
190
191
            labelList& snappedPoint
        ) const;

        //- Generate single point by interpolation or snapping
mattijs's avatar
mattijs committed
192
193
        template<class Type>
        Type generatePoint
mattijs's avatar
mattijs committed
194
        (
mattijs's avatar
mattijs committed
195
            const scalar s0,
mattijs's avatar
mattijs committed
196
            const Type& p0,
mattijs's avatar
mattijs committed
197
198
            const bool hasSnap0,
            const Type& snapP0,
mattijs's avatar
mattijs committed
199

mattijs's avatar
mattijs committed
200
            const scalar s1,
mattijs's avatar
mattijs committed
201
            const Type& p1,
mattijs's avatar
mattijs committed
202
203
            const bool hasSnap1,
            const Type& snapP1
mattijs's avatar
mattijs committed
204
205
        ) const;

mattijs's avatar
mattijs committed
206
        template<class Type>
mattijs's avatar
mattijs committed
207
208
209
        void generateTriPoints
        (
            const scalar s0,
mattijs's avatar
mattijs committed
210
            const Type& p0,
mattijs's avatar
mattijs committed
211
212
            const bool hasSnap0,
            const Type& snapP0,
mattijs's avatar
mattijs committed
213

mattijs's avatar
mattijs committed
214
            const scalar s1,
mattijs's avatar
mattijs committed
215
            const Type& p1,
mattijs's avatar
mattijs committed
216
217
            const bool hasSnap1,
            const Type& snapP1,
mattijs's avatar
mattijs committed
218

mattijs's avatar
mattijs committed
219
            const scalar s2,
mattijs's avatar
mattijs committed
220
            const Type& p2,
mattijs's avatar
mattijs committed
221
222
            const bool hasSnap2,
            const Type& snapP2,
mattijs's avatar
mattijs committed
223

mattijs's avatar
mattijs committed
224
            const scalar s3,
mattijs's avatar
mattijs committed
225
            const Type& p3,
mattijs's avatar
mattijs committed
226
227
            const bool hasSnap3,
            const Type& snapP3,
mattijs's avatar
mattijs committed
228
229
230
231
232

            DynamicList<Type>& points
        ) const;

        template<class Type>
mattijs's avatar
mattijs committed
233
        label generateFaceTriPoints
mattijs's avatar
mattijs committed
234
235
236
237
238
239
240
241
242
243
244
245
246
247
        (
            const volScalarField& cVals,
            const scalarField& pVals,

            const GeometricField<Type, fvPatchField, volMesh>& cCoords,
            const Field<Type>& pCoords,

            const DynamicList<Type>& snappedPoints,
            const labelList& snappedCc,
            const labelList& snappedPoint,
            const label faceI,

            const scalar neiVal,
            const Type& neiPt,
mattijs's avatar
mattijs committed
248
249
            const bool hasNeiSnap,
            const Type& neiSnapPt,
mattijs's avatar
mattijs committed
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269

            DynamicList<Type>& triPoints,
            DynamicList<label>& triMeshCells
        ) const;

        template<class Type>
        void generateTriPoints
        (
            const volScalarField& cVals,
            const scalarField& pVals,

            const GeometricField<Type, fvPatchField, volMesh>& cCoords,
            const Field<Type>& pCoords,

            const DynamicList<Type>& snappedPoints,
            const labelList& snappedCc,
            const labelList& snappedPoint,

            DynamicList<Type>& triPoints,
            DynamicList<label>& triMeshCells
mattijs's avatar
mattijs committed
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
        ) const;

        triSurface stitchTriPoints
        (
            const bool checkDuplicates,
            const List<point>& triPoints,
            labelList& triPointReverseMap,  // unmerged to merged point
            labelList& triMap               // merged to unmerged triangle
        ) const;

        //- Check single triangle for (topological) validity
        static bool validTri(const triSurface&, const label);

        //- Determine edge-face addressing
        void calcAddressing
        (
            const triSurface& surf,
            List<FixedList<label, 3> >& faceEdges,
            labelList& edgeFace0,
            labelList& edgeFace1,
            Map<labelList>& edgeFacesRest
        ) const;

        //- Determine orientation
        static void walkOrientation
        (
            const triSurface& surf,
            const List<FixedList<label, 3> >& faceEdges,
            const labelList& edgeFace0,
            const labelList& edgeFace1,
            const label seedTriI,
            labelList& flipState
        );

        //- Orient surface
        static void orientSurface
        (
            triSurface&,
            const List<FixedList<label, 3> >& faceEdges,
            const labelList& edgeFace0,
            const labelList& edgeFace1,
            const Map<labelList>& edgeFacesRest
        );
mattijs's avatar
mattijs committed
313

mattijs's avatar
mattijs committed
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
        //- Is triangle (given by 3 edges) not fully connected?
        static bool danglingTriangle
        (
            const FixedList<label, 3>& fEdges,
            const labelList& edgeFace1
        );

        //- Mark all non-fully connected triangles
        static label markDanglingTriangles
        (
            const List<FixedList<label, 3> >& faceEdges,
            const labelList& edgeFace0,
            const labelList& edgeFace1,
            const Map<labelList>& edgeFacesRest,
            boolList& keepTriangles
        );
mattijs's avatar
mattijs committed
330

mattijs's avatar
mattijs committed
331
332
333
334
        static triSurface subsetMesh
        (
            const triSurface& s,
            const labelList& newToOldFaces,
mattijs's avatar
mattijs committed
335
            labelList& oldToNewPoints,
mattijs's avatar
mattijs committed
336
            labelList& newToOldPoints
mattijs's avatar
mattijs committed
337
338
339
340
341
342
343
344
345
346
        );

public:

    //- Runtime type information
    TypeName("isoSurface");


    // Constructors

mattijs's avatar
mattijs committed
347
348
349
        //- Construct from cell values and point values. Uses boundaryField
        //  for boundary values. Requires on coupled patchfields to be set
        //  to the opposite cell value.
mattijs's avatar
mattijs committed
350
351
        isoSurface
        (
mattijs's avatar
mattijs committed
352
353
            const volScalarField& cellIsoVals,
            const scalarField& pointIsoVals,
mattijs's avatar
mattijs committed
354
355
356
            const scalar iso,
            const bool regularise,
            const scalar mergeTol = 1E-6    // fraction of bounding box
mattijs's avatar
mattijs committed
357
358
359
360
361
362
363
364
365
366
367
368
        );


    // Member Functions

        //- For every face original cell in mesh
        const labelList& meshCells() const
        {
            return meshCells_;
        }

        //- For every unmerged triangle point the point in the triSurface
mattijs's avatar
mattijs committed
369
        const labelList& triPointMergeMap() const
mattijs's avatar
mattijs committed
370
371
372
        {
            return triPointMergeMap_;
        }
mattijs's avatar
mattijs committed
373
374
375
376
377
378
379
380
381
382
383
384

        //- Interpolates cCoords,pCoords. Takes the original fields
        //  used to create the iso surface.
        template <class Type>
        tmp<Field<Type> > interpolate
        (
            const volScalarField& cellIsoVals,
            const scalarField& pointIsoVals,
            const GeometricField<Type, fvPatchField, volMesh>& cCoords,
            const Field<Type>& pCoords
        ) const;

mattijs's avatar
mattijs committed
385
386
387
388
389
390
391
392
393
};


// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

} // End namespace Foam

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

mattijs's avatar
mattijs committed
394
395
396
397
398
399
#ifdef NoRepository
#   include "isoSurfaceTemplates.C"
#endif

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

mattijs's avatar
mattijs committed
400
401
402
#endif

// ************************************************************************* //