isoSurface.H 14 KB
Newer Older
mattijs's avatar
mattijs committed
1 2 3 4
/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
OpenFOAM bot's avatar
OpenFOAM bot committed
5
    \\  /    A nd           | www.openfoam.com
OpenFOAM bot's avatar
OpenFOAM bot committed
6 7
     \\/     M anipulation  |
-------------------------------------------------------------------------------
OpenFOAM bot's avatar
OpenFOAM bot committed
8
    Copyright (C) 2011-2016 OpenFOAM Foundation
9
    Copyright (C) 2016-2019 OpenCFD Ltd.
mattijs's avatar
mattijs committed
10 11 12 13
-------------------------------------------------------------------------------
License
    This file is part of OpenFOAM.

14 15 16 17
    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 3 of the License, or
    (at your option) any later version.
mattijs's avatar
mattijs committed
18 19 20 21 22 23 24

    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
25
    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
mattijs's avatar
mattijs committed
26 27 28 29 30 31

Class
    Foam::isoSurface

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

    Note:
    - 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
39 40
    - regularisation does not do boundary points so as to preserve the
      boundary perfectly.
mattijs's avatar
mattijs committed
41 42 43
    - 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
44
    - on empty patches behaves like zero gradient.
mattijs's avatar
mattijs committed
45
    - does not do 2D correctly, creates non-flat iso surface.
mattijs's avatar
mattijs committed
46 47 48 49 50 51 52 53 54 55 56 57 58 59
    - on processor boundaries might two overlapping (identical) triangles
      (one from either side)

    The handling on coupled patches is a bit complex. All fields
    (values and coordinates) get rewritten so
    - empty patches get zerogradient (value) and facecentre (coordinates)
    - separated processor patch faces get interpolate (value) and facecentre
      (coordinates). (this is already the default for cyclics)
    - non-separated processor patch faces keep opposite value and cell centre

    Now the triangle generation on non-separated processor patch faces
    can use the neighbouring value. Any separated processor face or cyclic
    face gets handled just like any boundary face.

mattijs's avatar
mattijs committed
60 61
SourceFiles
    isoSurface.C
62
    isoSurfaceTemplates.C
mattijs's avatar
mattijs committed
63 64 65 66 67 68

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

#ifndef isoSurface_H
#define isoSurface_H

69
#include "bitSet.H"
mattijs's avatar
mattijs committed
70
#include "volFields.H"
mattijs's avatar
mattijs committed
71
#include "slicedVolFields.H"
72
#include "isoSurfaceBase.H"
mattijs's avatar
mattijs committed
73 74 75 76 77 78

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

namespace Foam
{

79 80
// Forward declarations

mattijs's avatar
mattijs committed
81
class fvMesh;
82 83
class plane;
class treeBoundBox;
84
class triSurface;
mattijs's avatar
mattijs committed
85 86 87 88 89 90 91

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

class isoSurface
:
92
    public isoSurfaceBase
mattijs's avatar
mattijs committed
93 94 95
{
    // Private data

mattijs's avatar
mattijs committed
96 97 98 99 100 101 102 103 104 105 106 107 108
        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
109 110


mattijs's avatar
mattijs committed
111
        //- Reference to mesh
mattijs's avatar
mattijs committed
112
        const fvMesh& mesh_;
mattijs's avatar
mattijs committed
113

mattijs's avatar
mattijs committed
114 115 116 117 118
        const scalarField& pVals_;

        //- Input volScalarField with separated coupled patches rewritten
        autoPtr<slicedVolScalarField> cValsPtr_;

mattijs's avatar
mattijs committed
119
        //- Regularise?
Mark Olesen's avatar
Mark Olesen committed
120
        const bool regularise_;
mattijs's avatar
mattijs committed
121

mattijs's avatar
mattijs committed
122 123 124
        //- When to merge points
        const scalar mergeDistance_;

mattijs's avatar
mattijs committed
125 126 127 128 129 130 131 132 133
        //- 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
134 135 136
        //- For every unmerged triangle point the point in the triSurface
        labelList triPointMergeMap_;

137 138 139 140
        //- triSurface points that have weighted interpolation
        DynamicList<label> interpolatedPoints_;

        //- corresponding original, unmerged points
141
        DynamicList<FixedList<label, 3>> interpolatedOldPoints_;
142 143

        //- corresponding weights
144
        DynamicList<FixedList<scalar, 3>> interpolationWeights_;
145

mattijs's avatar
mattijs committed
146 147 148

    // Private Member Functions

mattijs's avatar
mattijs committed
149 150 151 152 153 154
        // Point synchronisation

            //- Does tensor differ (to within mergeTolerance) from identity
            bool noTransform(const tensor& tt) const;

            //- Is patch a collocated (i.e. non-separated) coupled patch?
155
            static bool collocatedPatch(const polyPatch& pp);
mattijs's avatar
mattijs committed
156 157

            //- Per face whether is collocated
158
            bitSet collocatedFaces(const coupledPolyPatch&) const;
mattijs's avatar
mattijs committed
159

160
            //- Synchronise points on all non-separated coupled patches
mattijs's avatar
mattijs committed
161 162 163 164 165 166 167
            void syncUnseparatedPoints
            (
                pointField& collapsedPoint,
                const point& nullValue
            ) const;


mattijs's avatar
mattijs committed
168 169
        //- Get location of iso value as fraction inbetween s0,s1
        scalar isoFraction
mattijs's avatar
mattijs committed
170 171 172
        (
            const scalar s0,
            const scalar s1
mattijs's avatar
mattijs committed
173 174
        ) const;

mattijs's avatar
mattijs committed
175 176 177 178 179 180 181 182 183
        //- Check if any edge of a face is cut
        bool isEdgeOfFaceCut
        (
            const scalarField& pVals,
            const face& f,
            const bool ownLower,
            const bool neiLower
        ) const;

184
        //- Get neighbour value and position.
mattijs's avatar
mattijs committed
185 186 187
        void getNeighbour
        (
            const labelList& boundaryRegion,
mattijs's avatar
mattijs committed
188
            const volVectorField& meshC,
mattijs's avatar
mattijs committed
189
            const volScalarField& cVals,
190 191
            const label celli,
            const label facei,
mattijs's avatar
mattijs committed
192 193 194 195
            scalar& nbrValue,
            point& nbrPoint
        ) const;

196 197
        //- Determine for every face/cell whether it (possibly) generates
        //  triangles.
mattijs's avatar
mattijs committed
198
        void calcCutTypes
mattijs's avatar
mattijs committed
199
        (
mattijs's avatar
mattijs committed
200
            const labelList& boundaryRegion,
mattijs's avatar
mattijs committed
201
            const volVectorField& meshC,
mattijs's avatar
mattijs committed
202 203 204
            const volScalarField& cVals,
            const scalarField& pVals
        );
mattijs's avatar
mattijs committed
205 206 207 208 209 210

        static point calcCentre(const triSurface&);

        //- Determine per cc whether all near cuts can be snapped to single
        //  point.
        void calcSnappedCc
mattijs's avatar
mattijs committed
211
        (
mattijs's avatar
mattijs committed
212
            const labelList& boundaryRegion,
mattijs's avatar
mattijs committed
213
            const volVectorField& meshC,
mattijs's avatar
mattijs committed
214
            const volScalarField& cVals,
mattijs's avatar
mattijs committed
215
            const scalarField& pVals,
mattijs's avatar
mattijs committed
216
            DynamicList<point>& snappedPoints,
mattijs's avatar
mattijs committed
217 218 219 220 221 222 223
            labelList& snappedCc
        ) const;

        //- Determine per point whether all near cuts can be snapped to single
        //  point.
        void calcSnappedPoint
        (
224
            const bitSet& isBoundaryPoint,
mattijs's avatar
mattijs committed
225
            const labelList& boundaryRegion,
mattijs's avatar
mattijs committed
226
            const volVectorField& meshC,
mattijs's avatar
mattijs committed
227
            const volScalarField& cVals,
mattijs's avatar
mattijs committed
228
            const scalarField& pVals,
mattijs's avatar
mattijs committed
229
            DynamicList<point>& snappedPoints,
mattijs's avatar
mattijs committed
230 231 232
            labelList& snappedPoint
        ) const;

mattijs's avatar
mattijs committed
233 234 235 236

        //- Return input field with coupled (and empty) patch values rewritten
        template<class Type>
        tmp<SlicedGeometricField
237
        <Type, fvPatchField, slicedFvPatchField, volMesh>>
mattijs's avatar
mattijs committed
238 239 240 241 242
        adaptPatchFields
        (
            const GeometricField<Type, fvPatchField, volMesh>& fld
        ) const;

mattijs's avatar
mattijs committed
243
        //- Generate single point by interpolation or snapping
mattijs's avatar
mattijs committed
244 245
        template<class Type>
        Type generatePoint
mattijs's avatar
mattijs committed
246
        (
mattijs's avatar
mattijs committed
247
            const scalar s0,
mattijs's avatar
mattijs committed
248
            const Type& p0,
mattijs's avatar
mattijs committed
249 250
            const bool hasSnap0,
            const Type& snapP0,
mattijs's avatar
mattijs committed
251

mattijs's avatar
mattijs committed
252
            const scalar s1,
mattijs's avatar
mattijs committed
253
            const Type& p1,
mattijs's avatar
mattijs committed
254 255
            const bool hasSnap1,
            const Type& snapP1
mattijs's avatar
mattijs committed
256 257
        ) const;

258 259 260

        //- Note: cannot use simpler isoSurfaceCell::generateTriPoints since
        //  the need here to sometimes pass in remote 'snappedPoints'
mattijs's avatar
mattijs committed
261
        template<class Type>
mattijs's avatar
mattijs committed
262 263 264
        void generateTriPoints
        (
            const scalar s0,
mattijs's avatar
mattijs committed
265
            const Type& p0,
mattijs's avatar
mattijs committed
266 267
            const bool hasSnap0,
            const Type& snapP0,
mattijs's avatar
mattijs committed
268

mattijs's avatar
mattijs committed
269
            const scalar s1,
mattijs's avatar
mattijs committed
270
            const Type& p1,
mattijs's avatar
mattijs committed
271 272
            const bool hasSnap1,
            const Type& snapP1,
mattijs's avatar
mattijs committed
273

mattijs's avatar
mattijs committed
274
            const scalar s2,
mattijs's avatar
mattijs committed
275
            const Type& p2,
mattijs's avatar
mattijs committed
276 277
            const bool hasSnap2,
            const Type& snapP2,
mattijs's avatar
mattijs committed
278

mattijs's avatar
mattijs committed
279
            const scalar s3,
mattijs's avatar
mattijs committed
280
            const Type& p3,
mattijs's avatar
mattijs committed
281 282
            const bool hasSnap3,
            const Type& snapP3,
mattijs's avatar
mattijs committed
283

284
            DynamicList<Type>& pts
mattijs's avatar
mattijs committed
285 286 287
        ) const;

        template<class Type>
mattijs's avatar
mattijs committed
288
        label generateFaceTriPoints
mattijs's avatar
mattijs committed
289 290 291 292 293 294 295 296 297 298
        (
            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,
299
            const label facei,
mattijs's avatar
mattijs committed
300 301 302

            const scalar neiVal,
            const Type& neiPt,
mattijs's avatar
mattijs committed
303 304
            const bool hasNeiSnap,
            const Type& neiSnapPt,
mattijs's avatar
mattijs committed
305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324

            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
325 326
        ) const;

327
        template<class Type>
328
        static tmp<Field<Type>> interpolate
329 330 331 332
        (
            const label nPoints,
            const labelList& triPointMergeMap,
            const labelList& interpolatedPoints,
333 334
            const List<FixedList<label, 3>>& interpolatedOldPoints,
            const List<FixedList<scalar, 3>>& interpolationWeights,
335 336 337
            const DynamicList<Type>& unmergedValues
        );

mattijs's avatar
mattijs committed
338 339 340 341 342 343 344 345
        triSurface stitchTriPoints
        (
            const bool checkDuplicates,
            const List<point>& triPoints,
            labelList& triPointReverseMap,  // unmerged to merged point
            labelList& triMap               // merged to unmerged triangle
        ) const;

346 347
        //- Trim triangle to planes
        static void trimToPlanes
mattijs's avatar
mattijs committed
348
        (
349 350 351
            const PtrList<plane>& planes,
            const triPointRef& tri,
            DynamicList<point>& newTriPoints
mattijs's avatar
mattijs committed
352 353
        );

354 355
        //- Trim all triangles to box
        static void trimToBox
mattijs's avatar
mattijs committed
356
        (
357 358 359
            const treeBoundBox& bb,
            DynamicList<point>& triPoints,
            DynamicList<label>& triMeshCells
mattijs's avatar
mattijs committed
360
        );
mattijs's avatar
mattijs committed
361

362 363 364
        //- Trim all triangles to box. Determine interpolation
        //  for existing and new points
        static void trimToBox
mattijs's avatar
mattijs committed
365
        (
366 367 368 369 370
            const treeBoundBox& bb,
            DynamicList<point>& triPoints,
            DynamicList<label>& triMap,
            labelList& triPointMap,
            labelList& interpolatedPoints,
371 372
            List<FixedList<label, 3>>& interpolatedOldPoints,
            List<FixedList<scalar, 3>>& interpolationWeights
mattijs's avatar
mattijs committed
373 374 375 376
        );

        static triSurface subsetMesh
        (
377
            const triSurface&,
mattijs's avatar
mattijs committed
378
            const labelList& newToOldFaces,
mattijs's avatar
mattijs committed
379
            labelList& oldToNewPoints,
mattijs's avatar
mattijs committed
380
            labelList& newToOldPoints
mattijs's avatar
mattijs committed
381 382 383 384
        );

public:

385
    //- Declare friendship to share some functionality
386
    friend class isoSurfaceCell;
387
    friend class isoSurfaceTopo;
388

389 390 391
    //- Filtering type
    using isoSurfaceBase::filterType;

392

mattijs's avatar
mattijs committed
393 394 395 396 397 398
    //- Runtime type information
    TypeName("isoSurface");


    // Constructors

399 400 401 402 403 404 405 406 407 408 409 410 411 412 413 414 415
        //- Construct from cell values and point values.
        //  Uses boundaryField for boundary values.
        //  Holds reference to cellIsoVals and pointIsoVals.
        //
        //  \param bounds optional bounding box for trimming
        //  \param mergeTol fraction of mesh bounding box for merging points
        isoSurface
        (
            const volScalarField& cellValues,
            const scalarField& pointValues,
            const scalar iso,
            const isoSurfaceBase::filterType filter,
            const boundBox& bounds = boundBox::invertedBox,
            const scalar mergeTol = 1e-6
        );


416 417 418
        //- Construct from cell values and point values.
        //  Uses boundaryField for boundary values.
        //  Holds reference to cellIsoVals and pointIsoVals.
419 420 421
        //
        //  \param bounds optional bounding box for trimming
        //  \param mergeTol fraction of mesh bounding box for merging points
mattijs's avatar
mattijs committed
422 423
        isoSurface
        (
424 425
            const volScalarField& cellValues,
            const scalarField& pointValues,
mattijs's avatar
mattijs committed
426 427
            const scalar iso,
            const bool regularise,
428
            const boundBox& bounds = boundBox::invertedBox,
429
            const scalar mergeTol = 1e-6
mattijs's avatar
mattijs committed
430 431 432 433 434
        );


    // Member Functions

435 436 437
        //- Interpolates cCoords, pCoords.
        //  Uses the references to the original fields used to create the
        //  iso surface.
Henry's avatar
Henry committed
438
        template<class Type>
439
        tmp<Field<Type>> interpolate
mattijs's avatar
mattijs committed
440 441 442 443 444
        (
            const GeometricField<Type, fvPatchField, volMesh>& cCoords,
            const Field<Type>& pCoords
        ) const;

mattijs's avatar
mattijs committed
445 446 447 448 449 450 451 452 453
};


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

} // End namespace Foam

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

mattijs's avatar
mattijs committed
454
#ifdef NoRepository
455
    #include "isoSurfaceTemplates.C"
mattijs's avatar
mattijs committed
456 457 458 459
#endif

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

mattijs's avatar
mattijs committed
460 461 462
#endif

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