isoSurface.C 54.9 KB
Newer Older
mattijs's avatar
mattijs committed
1
2
3
4
/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
5
    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
mattijs's avatar
mattijs committed
6
7
8
9
10
     \\/     M anipulation  |
-------------------------------------------------------------------------------
License
    This file is part of OpenFOAM.

11
12
13
14
    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
15
16
17
18
19
20
21

    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
22
    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
mattijs's avatar
mattijs committed
23
24
25
26

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

#include "isoSurface.H"
mattijs's avatar
mattijs committed
27
#include "fvMesh.H"
mattijs's avatar
mattijs committed
28
29
#include "mergePoints.H"
#include "addToRunTimeSelectionTable.H"
mattijs's avatar
mattijs committed
30
#include "slicedVolFields.H"
mattijs's avatar
mattijs committed
31
#include "volFields.H"
mattijs's avatar
mattijs committed
32
#include "surfaceFields.H"
mattijs's avatar
mattijs committed
33
34
#include "OFstream.H"
#include "meshTools.H"
mattijs's avatar
mattijs committed
35
36
37

// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //

38
39
40
41
namespace Foam
{
defineTypeNameAndDebug(isoSurface, 0);
}
42

mattijs's avatar
mattijs committed
43
44
45

// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //

mattijs's avatar
mattijs committed
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
bool Foam::isoSurface::noTransform(const tensor& tt) const
{
    return
        (mag(tt.xx()-1) < mergeDistance_)
     && (mag(tt.yy()-1) < mergeDistance_)
     && (mag(tt.zz()-1) < mergeDistance_)
     && (mag(tt.xy()) < mergeDistance_)
     && (mag(tt.xz()) < mergeDistance_)
     && (mag(tt.yx()) < mergeDistance_)
     && (mag(tt.yz()) < mergeDistance_)
     && (mag(tt.zx()) < mergeDistance_)
     && (mag(tt.zy()) < mergeDistance_);
}


61
// Calculates per face whether couple is collocated.
mattijs's avatar
mattijs committed
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
bool Foam::isoSurface::collocatedPatch(const polyPatch& pp)
{
    const coupledPolyPatch& cpp = refCast<const coupledPolyPatch>(pp);
    return cpp.parallel() && !cpp.separated();
}


// Calculates per face whether couple is collocated.
Foam::PackedBoolList Foam::isoSurface::collocatedFaces
(
    const coupledPolyPatch& pp
) const
{
    // Initialise to false
    PackedBoolList collocated(pp.size());

78
    if (isA<processorPolyPatch>(pp))
mattijs's avatar
mattijs committed
79
    {
80
        if (collocatedPatch(pp))
mattijs's avatar
mattijs committed
81
        {
82
83
84
85
            forAll(pp, i)
            {
                collocated[i] = 1u;
            }
mattijs's avatar
mattijs committed
86
87
        }
    }
88
    else if (isA<cyclicPolyPatch>(pp))
mattijs's avatar
mattijs committed
89
    {
90
        if (collocatedPatch(pp))
mattijs's avatar
mattijs committed
91
        {
92
93
94
95
            forAll(pp, i)
            {
                collocated[i] = 1u;
            }
mattijs's avatar
mattijs committed
96
97
        }
    }
98
    else
mattijs's avatar
mattijs committed
99
    {
100
101
102
103
104
        FatalErrorIn
        (
            "isoSurface::collocatedFaces(const coupledPolyPatch&) const"
        )   << "Unhandled coupledPolyPatch type " << pp.type()
            << abort(FatalError);
mattijs's avatar
mattijs committed
105
    }
106
    return collocated;
mattijs's avatar
mattijs committed
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
}


void Foam::isoSurface::syncUnseparatedPoints
(
    pointField& pointValues,
    const point& nullValue
) const
{
    // Until syncPointList handles separated coupled patches with multiple
    // transforms do our own synchronisation of non-separated patches only
    const polyBoundaryMesh& patches = mesh_.boundaryMesh();

    if (Pstream::parRun())
    {
        // Send
        forAll(patches, patchI)
        {
            if
            (
                isA<processorPolyPatch>(patches[patchI])
             && patches[patchI].nPoints() > 0
mattijs's avatar
mattijs committed
129
             && collocatedPatch(patches[patchI])
mattijs's avatar
mattijs committed
130
131
132
133
134
            )
            {
                const processorPolyPatch& pp =
                    refCast<const processorPolyPatch>(patches[patchI]);

mattijs's avatar
mattijs committed
135
136
137
138
                const labelList& meshPts = pp.meshPoints();
                const labelList& nbrPts = pp.neighbPoints();

                pointField patchInfo(meshPts.size());
mattijs's avatar
mattijs committed
139

mattijs's avatar
mattijs committed
140
                forAll(nbrPts, pointI)
mattijs's avatar
mattijs committed
141
                {
mattijs's avatar
mattijs committed
142
143
                    label nbrPointI = nbrPts[pointI];
                    patchInfo[nbrPointI] = pointValues[meshPts[pointI]];
mattijs's avatar
mattijs committed
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
                }

                OPstream toNbr(Pstream::blocking, pp.neighbProcNo());
                toNbr << patchInfo;
            }
        }

        // Receive and combine.

        forAll(patches, patchI)
        {
            if
            (
                isA<processorPolyPatch>(patches[patchI])
             && patches[patchI].nPoints() > 0
mattijs's avatar
mattijs committed
159
             && collocatedPatch(patches[patchI])
mattijs's avatar
mattijs committed
160
161
162
163
164
            )
            {
                const processorPolyPatch& pp =
                    refCast<const processorPolyPatch>(patches[patchI]);

mattijs's avatar
mattijs committed
165
                pointField nbrPatchInfo(pp.nPoints());
mattijs's avatar
mattijs committed
166
167
168
169
170
171
172
                {
                    // We do not know the number of points on the other side
                    // so cannot use Pstream::read.
                    IPstream fromNbr(Pstream::blocking, pp.neighbProcNo());
                    fromNbr >> nbrPatchInfo;
                }

mattijs's avatar
mattijs committed
173
174
175
                const labelList& meshPts = pp.meshPoints();

                forAll(meshPts, pointI)
176
                {
mattijs's avatar
mattijs committed
177
178
179
180
181
182
                    label meshPointI = meshPts[pointI];
                    minEqOp<point>()
                    (
                        pointValues[meshPointI],
                        nbrPatchInfo[pointI]
                    );
mattijs's avatar
mattijs committed
183
184
185
186
187
                }
            }
        }
    }

188
189
190
191
192
193
194
    // Do the cyclics.
    forAll(patches, patchI)
    {
        if (isA<cyclicPolyPatch>(patches[patchI]))
        {
            const cyclicPolyPatch& cycPatch =
                refCast<const cyclicPolyPatch>(patches[patchI]);
mattijs's avatar
mattijs committed
195

196
197
198
            if (cycPatch.owner() && collocatedPatch(cycPatch))
            {
                // Owner does all.
mattijs's avatar
mattijs committed
199

200
201
202
203
                const edgeList& coupledPoints = cycPatch.coupledPoints();
                const labelList& meshPts = cycPatch.meshPoints();
                const cyclicPolyPatch& nbrPatch = cycPatch.neighbPatch();
                const labelList& nbrMeshPoints = nbrPatch.meshPoints();
mattijs's avatar
mattijs committed
204

205
206
                pointField half0Values(coupledPoints.size());
                pointField half1Values(coupledPoints.size());
mattijs's avatar
mattijs committed
207

208
209
210
211
212
213
214
215
216
217
218
219
                forAll(coupledPoints, i)
                {
                    const edge& e = coupledPoints[i];
                    half0Values[i] = pointValues[meshPts[e[0]]];
                    half1Values[i] = pointValues[nbrMeshPoints[e[1]]];
                }

                forAll(coupledPoints, i)
                {
                    const edge& e = coupledPoints[i];
                    label p0 = meshPts[e[0]];
                    label p1 = nbrMeshPoints[e[1]];
mattijs's avatar
mattijs committed
220

221
222
223
224
225
226
227
228
229
230
231
                    minEqOp<point>()(pointValues[p0], half1Values[i]);
                    minEqOp<point>()(pointValues[p1], half0Values[i]);
                }
            }
        }
    }

    // Synchronize multiple shared points.
    const globalMeshData& pd = mesh_.globalData();

    if (pd.nGlobalPoints() > 0)
mattijs's avatar
mattijs committed
232
    {
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
        // Values on shared points.
        pointField sharedPts(pd.nGlobalPoints(), nullValue);

        forAll(pd.sharedPointLabels(), i)
        {
            label meshPointI = pd.sharedPointLabels()[i];
            // Fill my entries in the shared points
            sharedPts[pd.sharedPointAddr()[i]] = pointValues[meshPointI];
        }

        // Combine on master.
        Pstream::listCombineGather(sharedPts, minEqOp<point>());
        Pstream::listCombineScatter(sharedPts);

        // Now we will all have the same information. Merge it back with
        // my local information.
        forAll(pd.sharedPointLabels(), i)
        {
            label meshPointI = pd.sharedPointLabels()[i];
            pointValues[meshPointI] = sharedPts[pd.sharedPointAddr()[i]];
        }
mattijs's avatar
mattijs committed
254
255
256
257
    }
}


mattijs's avatar
mattijs committed
258
259
260
261
262
263
264
Foam::scalar Foam::isoSurface::isoFraction
(
    const scalar s0,
    const scalar s1
) const
{
    scalar d = s1-s0;
mattijs's avatar
mattijs committed
265

mattijs's avatar
mattijs committed
266
267
268
269
270
271
272
273
274
    if (mag(d) > VSMALL)
    {
        return (iso_-s0)/d;
    }
    else
    {
        return -1.0;
    }
}
mattijs's avatar
mattijs committed
275

mattijs's avatar
mattijs committed
276

mattijs's avatar
mattijs committed
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
bool Foam::isoSurface::isEdgeOfFaceCut
(
    const scalarField& pVals,
    const face& f,
    const bool ownLower,
    const bool neiLower
) const
{
    forAll(f, fp)
    {
        bool fpLower = (pVals[f[fp]] < iso_);
        if
        (
            (fpLower != ownLower)
         || (fpLower != neiLower)
         || (fpLower != (pVals[f[f.fcIndex(fp)]] < iso_))
        )
        {
            return true;
        }
    }
    return false;
}


mattijs's avatar
mattijs committed
302
303
304
305
// Get neighbour value and position.
void Foam::isoSurface::getNeighbour
(
    const labelList& boundaryRegion,
mattijs's avatar
mattijs committed
306
    const volVectorField& meshC,
mattijs's avatar
mattijs committed
307
308
309
310
311
312
313
314
315
316
317
318
319
320
    const volScalarField& cVals,
    const label cellI,
    const label faceI,
    scalar& nbrValue,
    point& nbrPoint
) const
{
    const labelList& own = mesh_.faceOwner();
    const labelList& nei = mesh_.faceNeighbour();

    if (mesh_.isInternalFace(faceI))
    {
        label nbr = (own[faceI] == cellI ? nei[faceI] : own[faceI]);
        nbrValue = cVals[nbr];
mattijs's avatar
mattijs committed
321
        nbrPoint = meshC[nbr];
mattijs's avatar
mattijs committed
322
323
324
325
326
327
328
329
    }
    else
    {
        label bFaceI = faceI-mesh_.nInternalFaces();
        label patchI = boundaryRegion[bFaceI];
        const polyPatch& pp = mesh_.boundaryMesh()[patchI];
        label patchFaceI = faceI-pp.start();

mattijs's avatar
mattijs committed
330
331
        nbrValue = cVals.boundaryField()[patchI][patchFaceI];
        nbrPoint = meshC.boundaryField()[patchI][patchFaceI];
mattijs's avatar
mattijs committed
332
333
334
335
    }
}


mattijs's avatar
mattijs committed
336
337
// Determine for every face/cell whether it (possibly) generates triangles.
void Foam::isoSurface::calcCutTypes
mattijs's avatar
mattijs committed
338
(
mattijs's avatar
mattijs committed
339
    const labelList& boundaryRegion,
mattijs's avatar
mattijs committed
340
    const volVectorField& meshC,
mattijs's avatar
mattijs committed
341
342
343
    const volScalarField& cVals,
    const scalarField& pVals
)
mattijs's avatar
mattijs committed
344
{
mattijs's avatar
mattijs committed
345
346
347
348
349
350
    const polyBoundaryMesh& patches = mesh_.boundaryMesh();
    const labelList& own = mesh_.faceOwner();
    const labelList& nei = mesh_.faceNeighbour();

    faceCutType_.setSize(mesh_.nFaces());
    faceCutType_ = NOTCUT;
mattijs's avatar
mattijs committed
351

mattijs's avatar
mattijs committed
352
    for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
mattijs's avatar
mattijs committed
353
    {
mattijs's avatar
mattijs committed
354
355
        // CC edge.
        bool ownLower = (cVals[own[faceI]] < iso_);
mattijs's avatar
mattijs committed
356
357
358
359
360
361

        scalar nbrValue;
        point nbrPoint;
        getNeighbour
        (
            boundaryRegion,
mattijs's avatar
mattijs committed
362
            meshC,
mattijs's avatar
mattijs committed
363
364
365
366
367
368
369
370
            cVals,
            own[faceI],
            faceI,
            nbrValue,
            nbrPoint
        );

        bool neiLower = (nbrValue < iso_);
mattijs's avatar
mattijs committed
371
372

        if (ownLower != neiLower)
mattijs's avatar
mattijs committed
373
        {
mattijs's avatar
mattijs committed
374
375
376
377
            faceCutType_[faceI] = CUT;
        }
        else
        {
mattijs's avatar
mattijs committed
378
379
            // See if any mesh edge is cut by looping over all the edges of the
            // face.
mattijs's avatar
mattijs committed
380
            const face f = mesh_.faces()[faceI];
mattijs's avatar
mattijs committed
381

mattijs's avatar
mattijs committed
382
            if (isEdgeOfFaceCut(pVals, f, ownLower, neiLower))
mattijs's avatar
mattijs committed
383
            {
mattijs's avatar
mattijs committed
384
                faceCutType_[faceI] = CUT;
mattijs's avatar
mattijs committed
385
386
387
            }
        }
    }
mattijs's avatar
mattijs committed
388

mattijs's avatar
mattijs committed
389
    forAll(patches, patchI)
mattijs's avatar
mattijs committed
390
    {
mattijs's avatar
mattijs committed
391
        const polyPatch& pp = patches[patchI];
mattijs's avatar
mattijs committed
392

mattijs's avatar
mattijs committed
393
        label faceI = pp.start();
mattijs's avatar
mattijs committed
394

mattijs's avatar
mattijs committed
395
        forAll(pp, i)
mattijs's avatar
mattijs committed
396
        {
mattijs's avatar
mattijs committed
397
            bool ownLower = (cVals[own[faceI]] < iso_);
mattijs's avatar
mattijs committed
398

mattijs's avatar
mattijs committed
399
400
401
402
403
            scalar nbrValue;
            point nbrPoint;
            getNeighbour
            (
                boundaryRegion,
mattijs's avatar
mattijs committed
404
                meshC,
mattijs's avatar
mattijs committed
405
406
407
408
409
410
                cVals,
                own[faceI],
                faceI,
                nbrValue,
                nbrPoint
            );
mattijs's avatar
mattijs committed
411

mattijs's avatar
mattijs committed
412
            bool neiLower = (nbrValue < iso_);
mattijs's avatar
mattijs committed
413

mattijs's avatar
mattijs committed
414
415
416
417
418
419
420
            if (ownLower != neiLower)
            {
                faceCutType_[faceI] = CUT;
            }
            else
            {
                // Mesh edge.
mattijs's avatar
mattijs committed
421
                const face f = mesh_.faces()[faceI];
mattijs's avatar
mattijs committed
422

mattijs's avatar
mattijs committed
423
                if (isEdgeOfFaceCut(pVals, f, ownLower, neiLower))
mattijs's avatar
mattijs committed
424
                {
mattijs's avatar
mattijs committed
425
                    faceCutType_[faceI] = CUT;
mattijs's avatar
mattijs committed
426
                }
mattijs's avatar
mattijs committed
427
            }
mattijs's avatar
mattijs committed
428

mattijs's avatar
mattijs committed
429
            faceI++;
mattijs's avatar
mattijs committed
430
        }
mattijs's avatar
mattijs committed
431
    }
mattijs's avatar
mattijs committed
432
433
434



mattijs's avatar
mattijs committed
435
436
437
    nCutCells_ = 0;
    cellCutType_.setSize(mesh_.nCells());
    cellCutType_ = NOTCUT;
mattijs's avatar
mattijs committed
438

mattijs's avatar
mattijs committed
439
440
441
442
443
    for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
    {
        if (faceCutType_[faceI] != NOTCUT)
        {
            if (cellCutType_[own[faceI]] == NOTCUT)
mattijs's avatar
mattijs committed
444
            {
mattijs's avatar
mattijs committed
445
446
                cellCutType_[own[faceI]] = CUT;
                nCutCells_++;
mattijs's avatar
mattijs committed
447
            }
mattijs's avatar
mattijs committed
448
            if (cellCutType_[nei[faceI]] == NOTCUT)
mattijs's avatar
mattijs committed
449
            {
mattijs's avatar
mattijs committed
450
451
                cellCutType_[nei[faceI]] = CUT;
                nCutCells_++;
mattijs's avatar
mattijs committed
452
453
            }
        }
mattijs's avatar
mattijs committed
454
455
456
457
    }
    for (label faceI = mesh_.nInternalFaces(); faceI < mesh_.nFaces(); faceI++)
    {
        if (faceCutType_[faceI] != NOTCUT)
mattijs's avatar
mattijs committed
458
        {
mattijs's avatar
mattijs committed
459
460
461
462
463
            if (cellCutType_[own[faceI]] == NOTCUT)
            {
                cellCutType_[own[faceI]] = CUT;
                nCutCells_++;
            }
mattijs's avatar
mattijs committed
464
465
        }
    }
mattijs's avatar
mattijs committed
466
467
468
469

    if (debug)
    {
        Pout<< "isoSurface : detected " << nCutCells_
mattijs's avatar
mattijs committed
470
471
            << " candidate cut cells (out of " << mesh_.nCells()
            << ")." << endl;
mattijs's avatar
mattijs committed
472
    }
mattijs's avatar
mattijs committed
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
}


// Return the two common points between two triangles
Foam::labelPair Foam::isoSurface::findCommonPoints
(
    const labelledTri& tri0,
    const labelledTri& tri1
)
{
    labelPair common(-1, -1);

    label fp0 = 0;
    label fp1 = findIndex(tri1, tri0[fp0]);

    if (fp1 == -1)
    {
        fp0 = 1;
        fp1 = findIndex(tri1, tri0[fp0]);
    }

    if (fp1 != -1)
    {
        // So tri0[fp0] is tri1[fp1]

        // Find next common point
        label fp0p1 = tri0.fcIndex(fp0);
        label fp1p1 = tri1.fcIndex(fp1);
        label fp1m1 = tri1.rcIndex(fp1);

        if (tri0[fp0p1] == tri1[fp1p1] || tri0[fp0p1] == tri1[fp1m1])
        {
            common[0] = tri0[fp0];
            common[1] = tri0[fp0p1];
        }
    }
    return common;
}


// Caculate centre of surface.
Foam::point Foam::isoSurface::calcCentre(const triSurface& s)
{
    vector sum = vector::zero;

    forAll(s, i)
    {
        sum += s[i].centre(s.points());
    }
    return sum/s.size();
}


// Replace surface (localPoints, localTris) with single point. Returns
// point. Destructs arguments.
Foam::pointIndexHit Foam::isoSurface::collapseSurface
(
    pointField& localPoints,
    DynamicList<labelledTri, 64>& localTris
)
{
    pointIndexHit info(false, vector::zero, localTris.size());

536
    if (localTris.size() == 1)
mattijs's avatar
mattijs committed
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
    {
        const labelledTri& tri = localTris[0];
        info.setPoint(tri.centre(localPoints));
        info.setHit();
    }
    else if (localTris.size() == 2)
    {
        // Check if the two triangles share an edge.
        const labelledTri& tri0 = localTris[0];
        const labelledTri& tri1 = localTris[0];

        labelPair shared = findCommonPoints(tri0, tri1);

        if (shared[0] != -1)
        {
            info.setPoint
mattijs's avatar
mattijs committed
553
            (
mattijs's avatar
mattijs committed
554
555
556
557
558
                0.5
              * (
                    tri0.centre(localPoints)
                  + tri1.centre(localPoints)
                )
mattijs's avatar
mattijs committed
559
            );
mattijs's avatar
mattijs committed
560
            info.setHit();
mattijs's avatar
mattijs committed
561
        }
mattijs's avatar
mattijs committed
562
    }
563
    else if (localTris.size())
mattijs's avatar
mattijs committed
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
    {
        // Check if single region. Rare situation.
        triSurface surf
        (
            localTris,
            geometricSurfacePatchList(0),
            localPoints,
            true
        );
        localTris.clearStorage();

        labelList faceZone;
        label nZones = surf.markZones
        (
            boolList(surf.nEdges(), false),
            faceZone
        );

        if (nZones == 1)
        {
            info.setPoint(calcCentre(surf));
            info.setHit();
        }
    }

    return info;
}


mattijs's avatar
mattijs committed
593
594
// Determine per cell centre whether all the intersections get collapsed
// to a single point
mattijs's avatar
mattijs committed
595
596
void Foam::isoSurface::calcSnappedCc
(
mattijs's avatar
mattijs committed
597
    const labelList& boundaryRegion,
mattijs's avatar
mattijs committed
598
    const volVectorField& meshC,
mattijs's avatar
mattijs committed
599
    const volScalarField& cVals,
mattijs's avatar
mattijs committed
600
601
602
603
604
605
606
    const scalarField& pVals,

    DynamicList<point>& snappedPoints,
    labelList& snappedCc
) const
{
    const pointField& pts = mesh_.points();
mattijs's avatar
mattijs committed
607
    const pointField& cc = mesh_.cellCentres();
mattijs's avatar
mattijs committed
608
609
610
611
612

    snappedCc.setSize(mesh_.nCells());
    snappedCc = -1;

    // Work arrays
mattijs's avatar
mattijs committed
613
    DynamicList<point, 64> localTriPoints(64);
mattijs's avatar
mattijs committed
614
615
616

    forAll(mesh_.cells(), cellI)
    {
mattijs's avatar
mattijs committed
617
        if (cellCutType_[cellI] == CUT)
mattijs's avatar
mattijs committed
618
        {
mattijs's avatar
mattijs committed
619
620
621
622
            scalar cVal = cVals[cellI];

            const cell& cFaces = mesh_.cells()[cellI];

mattijs's avatar
mattijs committed
623
624
625
            localTriPoints.clear();
            label nOther = 0;
            point otherPointSum = vector::zero;
mattijs's avatar
mattijs committed
626
627
628

            // Create points for all intersections close to cell centre
            // (i.e. from pyramid edges)
mattijs's avatar
mattijs committed
629
630
631

            forAll(cFaces, cFaceI)
            {
mattijs's avatar
mattijs committed
632
633
634
635
636
637
638
                label faceI = cFaces[cFaceI];

                scalar nbrValue;
                point nbrPoint;
                getNeighbour
                (
                    boundaryRegion,
mattijs's avatar
mattijs committed
639
                    meshC,
mattijs's avatar
mattijs committed
640
641
642
643
644
645
646
647
648
649
650
651
652
                    cVals,
                    cellI,
                    faceI,
                    nbrValue,
                    nbrPoint
                );

                // Calculate intersection points of edges to cell centre
                FixedList<scalar, 3> s;
                FixedList<point, 3> pt;

                // From cc to neighbour cc.
                s[2] = isoFraction(cVal, nbrValue);
mattijs's avatar
mattijs committed
653
                pt[2] = (1.0-s[2])*cc[cellI] + s[2]*nbrPoint;
mattijs's avatar
mattijs committed
654

mattijs's avatar
mattijs committed
655
656
657
658
                const face& f = mesh_.faces()[cFaces[cFaceI]];

                forAll(f, fp)
                {
mattijs's avatar
mattijs committed
659
660
661
                    // From cc to fp
                    label p0 = f[fp];
                    s[0] = isoFraction(cVal, pVals[p0]);
mattijs's avatar
mattijs committed
662
                    pt[0] = (1.0-s[0])*cc[cellI] + s[0]*pts[p0];
mattijs's avatar
mattijs committed
663
664
665
666

                    // From cc to fp+1
                    label p1 = f[f.fcIndex(fp)];
                    s[1] = isoFraction(cVal, pVals[p1]);
mattijs's avatar
mattijs committed
667
                    pt[1] = (1.0-s[1])*cc[cellI] + s[1]*pts[p1];
mattijs's avatar
mattijs committed
668
669
670
671
672
673
674
675
676
677
678
679
680

                    if
                    (
                        (s[0] >= 0.0 && s[0] <= 0.5)
                     && (s[1] >= 0.0 && s[1] <= 0.5)
                     && (s[2] >= 0.0 && s[2] <= 0.5)
                    )
                    {
                        localTriPoints.append(pt[0]);
                        localTriPoints.append(pt[1]);
                        localTriPoints.append(pt[2]);
                    }
                    else
mattijs's avatar
mattijs committed
681
                    {
mattijs's avatar
mattijs committed
682
683
                        // Get average of all other points
                        forAll(s, i)
mattijs's avatar
mattijs committed
684
                        {
mattijs's avatar
mattijs committed
685
686
687
688
689
                            if (s[i] >= 0.0 && s[i] <= 0.5)
                            {
                                otherPointSum += pt[i];
                                nOther++;
                            }
mattijs's avatar
mattijs committed
690
691
692
693
694
                        }
                    }
                }
            }

mattijs's avatar
mattijs committed
695
            if (localTriPoints.size() == 0)
mattijs's avatar
mattijs committed
696
            {
mattijs's avatar
mattijs committed
697
698
699
700
701
702
                // No complete triangles. Get average of all intersection
                // points.
                if (nOther > 0)
                {
                    snappedCc[cellI] = snappedPoints.size();
                    snappedPoints.append(otherPointSum/nOther);
mattijs's avatar
mattijs committed
703

mattijs's avatar
mattijs committed
704
705
706
707
                    //Pout<< "    point:" << pointI
                    //    << " replacing coord:" << mesh_.points()[pointI]
                    //    << " by average:" << collapsedPoint[pointI] << endl;
                }
mattijs's avatar
mattijs committed
708
            }
mattijs's avatar
mattijs committed
709
            else if (localTriPoints.size() == 3)
mattijs's avatar
mattijs committed
710
            {
mattijs's avatar
mattijs committed
711
712
713
                // Single triangle. No need for any analysis. Average points.
                pointField points;
                points.transfer(localTriPoints);
mattijs's avatar
mattijs committed
714
                snappedCc[cellI] = snappedPoints.size();
mattijs's avatar
mattijs committed
715
                snappedPoints.append(sum(points)/points.size());
mattijs's avatar
mattijs committed
716

mattijs's avatar
mattijs committed
717
718
719
                //Pout<< "    point:" << pointI
                //    << " replacing coord:" << mesh_.points()[pointI]
                //    << " by average:" << collapsedPoint[pointI] << endl;
mattijs's avatar
mattijs committed
720
721
722
            }
            else
            {
mattijs's avatar
mattijs committed
723
                // Convert points into triSurface.
mattijs's avatar
mattijs committed
724

mattijs's avatar
mattijs committed
725
726
727
728
729
730
731
732
733
                // Merge points and compact out non-valid triangles
                labelList triMap;               // merged to unmerged triangle
                labelList triPointReverseMap;   // unmerged to merged point
                triSurface surf
                (
                    stitchTriPoints
                    (
                        false,              // do not check for duplicate tris
                        localTriPoints,
734
                        triPointReverseMap,
mattijs's avatar
mattijs committed
735
736
737
                        triMap
                    )
                );
mattijs's avatar
mattijs committed
738

mattijs's avatar
mattijs committed
739
740
741
742
743
744
                labelList faceZone;
                label nZones = surf.markZones
                (
                    boolList(surf.nEdges(), false),
                    faceZone
                );
mattijs's avatar
mattijs committed
745

mattijs's avatar
mattijs committed
746
                if (nZones == 1)
mattijs's avatar
mattijs committed
747
                {
mattijs's avatar
mattijs committed
748
                    snappedCc[cellI] = snappedPoints.size();
mattijs's avatar
mattijs committed
749
750
751
752
                    snappedPoints.append(calcCentre(surf));
                    //Pout<< "    point:" << pointI << " nZones:" << nZones
                    //    << " replacing coord:" << mesh_.points()[pointI]
                    //    << " by average:" << collapsedPoint[pointI] << endl;
mattijs's avatar
mattijs committed
753
754
755
                }
            }
        }
mattijs's avatar
mattijs committed
756
757
758
759
    }
}


mattijs's avatar
mattijs committed
760
761
// Determine per meshpoint whether all the intersections get collapsed
// to a single point
mattijs's avatar
mattijs committed
762
763
void Foam::isoSurface::calcSnappedPoint
(
764
    const PackedBoolList& isBoundaryPoint,
mattijs's avatar
mattijs committed
765
    const labelList& boundaryRegion,
mattijs's avatar
mattijs committed
766
    const volVectorField& meshC,
mattijs's avatar
mattijs committed
767
    const volScalarField& cVals,
mattijs's avatar
mattijs committed
768
769
770
771
772
773
    const scalarField& pVals,

    DynamicList<point>& snappedPoints,
    labelList& snappedPoint
) const
{
mattijs's avatar
mattijs committed
774
    const pointField& pts = mesh_.points();
mattijs's avatar
mattijs committed
775
    const pointField& cc = mesh_.cellCentres();
mattijs's avatar
mattijs committed
776

777
    pointField collapsedPoint(mesh_.nPoints(), point::max);
mattijs's avatar
mattijs committed
778

mattijs's avatar
mattijs committed
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797

    // Work arrays
    DynamicList<point, 64> localTriPoints(100);

    forAll(mesh_.pointFaces(), pointI)
    {
        if (isBoundaryPoint.get(pointI) == 1)
        {
            continue;
        }

        const labelList& pFaces = mesh_.pointFaces()[pointI];

        bool anyCut = false;

        forAll(pFaces, i)
        {
            label faceI = pFaces[i];

mattijs's avatar
mattijs committed
798
            if (faceCutType_[faceI] == CUT)
mattijs's avatar
mattijs committed
799
800
801
802
803
804
805
806
807
808
809
810
811
            {
                anyCut = true;
                break;
            }
        }

        if (!anyCut)
        {
            continue;
        }


        localTriPoints.clear();
mattijs's avatar
mattijs committed
812
813
        label nOther = 0;
        point otherPointSum = vector::zero;
mattijs's avatar
mattijs committed
814
815
816

        forAll(pFaces, pFaceI)
        {
mattijs's avatar
mattijs committed
817
818
819
            // Create points for all intersections close to point
            // (i.e. from pyramid edges)

mattijs's avatar
mattijs committed
820
821
822
823
            label faceI = pFaces[pFaceI];
            const face& f = mesh_.faces()[faceI];
            label own = mesh_.faceOwner()[faceI];

mattijs's avatar
mattijs committed
824
            // Get neighbour value
mattijs's avatar
mattijs committed
825
826
827
828
829
            scalar nbrValue;
            point nbrPoint;
            getNeighbour
            (
                boundaryRegion,
mattijs's avatar
mattijs committed
830
                meshC,
mattijs's avatar
mattijs committed
831
832
833
834
835
836
837
838
839
840
841
842
843
                cVals,
                own,
                faceI,
                nbrValue,
                nbrPoint
            );

            // Calculate intersection points of edges emanating from point
            FixedList<scalar, 4> s;
            FixedList<point, 4> pt;

            label fp = findIndex(f, pointI);
            s[0] = isoFraction(pVals[pointI], cVals[own]);
mattijs's avatar
mattijs committed
844
            pt[0] = (1.0-s[0])*pts[pointI] + s[0]*cc[own];
mattijs's avatar
mattijs committed
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862

            s[1] = isoFraction(pVals[pointI], nbrValue);
            pt[1] = (1.0-s[1])*pts[pointI] + s[1]*nbrPoint;

            label nextPointI = f[f.fcIndex(fp)];
            s[2] = isoFraction(pVals[pointI], pVals[nextPointI]);
            pt[2] = (1.0-s[2])*pts[pointI] + s[2]*pts[nextPointI];

            label prevPointI = f[f.rcIndex(fp)];
            s[3] = isoFraction(pVals[pointI], pVals[prevPointI]);
            pt[3] = (1.0-s[3])*pts[pointI] + s[3]*pts[prevPointI];

            if
            (
                (s[0] >= 0.0 && s[0] <= 0.5)
             && (s[1] >= 0.0 && s[1] <= 0.5)
             && (s[2] >= 0.0 && s[2] <= 0.5)
            )
mattijs's avatar
mattijs committed
863
            {
mattijs's avatar
mattijs committed
864
865
866
                localTriPoints.append(pt[0]);
                localTriPoints.append(pt[1]);
                localTriPoints.append(pt[2]);
mattijs's avatar
mattijs committed
867
            }
mattijs's avatar
mattijs committed
868
869
870
871
872
873
            if
            (
                (s[0] >= 0.0 && s[0] <= 0.5)
             && (s[1] >= 0.0 && s[1] <= 0.5)
             && (s[3] >= 0.0 && s[3] <= 0.5)
            )
mattijs's avatar
mattijs committed
874
            {
mattijs's avatar
mattijs committed
875
876
877
                localTriPoints.append(pt[3]);
                localTriPoints.append(pt[0]);
                localTriPoints.append(pt[1]);
mattijs's avatar
mattijs committed
878
879
            }

mattijs's avatar
mattijs committed
880
881
            // Get average of points as fallback
            forAll(s, i)
mattijs's avatar
mattijs committed
882
            {
mattijs's avatar
mattijs committed
883
                if (s[i] >= 0.0 && s[i] <= 0.5)
mattijs's avatar
mattijs committed
884
                {
mattijs's avatar
mattijs committed
885
886
                    otherPointSum += pt[i];
                    nOther++;
mattijs's avatar
mattijs committed
887
888
889
890
891
892
                }
            }
        }

        if (localTriPoints.size() == 0)
        {
mattijs's avatar
mattijs committed
893
894
895
896
897
898
            // No complete triangles. Get average of all intersection
            // points.
            if (nOther > 0)
            {
                collapsedPoint[pointI] = otherPointSum/nOther;
            }
mattijs's avatar
mattijs committed
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
        }
        else if (localTriPoints.size() == 3)
        {
            // Single triangle. No need for any analysis. Average points.
            pointField points;
            points.transfer(localTriPoints);
            collapsedPoint[pointI] = sum(points)/points.size();
        }
        else
        {
            // Convert points into triSurface.

            // Merge points and compact out non-valid triangles
            labelList triMap;               // merged to unmerged triangle
            labelList triPointReverseMap;   // unmerged to merged point
            triSurface surf
            (
                stitchTriPoints
                (
                    false,                  // do not check for duplicate tris
                    localTriPoints,
920
                    triPointReverseMap,
mattijs's avatar
mattijs committed
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
                    triMap
                )
            );

            labelList faceZone;
            label nZones = surf.markZones
            (
                boolList(surf.nEdges(), false),
                faceZone
            );

            if (nZones == 1)
            {
                collapsedPoint[pointI] = calcCentre(surf);
            }
        }
    }

mattijs's avatar
mattijs committed
939

mattijs's avatar
mattijs committed
940
    // Synchronise snap location
941
    syncUnseparatedPoints(collapsedPoint, point::max);
mattijs's avatar
mattijs committed
942

mattijs's avatar
mattijs committed
943
944
945
946
947
948

    snappedPoint.setSize(mesh_.nPoints());
    snappedPoint = -1;

    forAll(collapsedPoint, pointI)
    {
949
        if (collapsedPoint[pointI] != point::max)
mattijs's avatar
mattijs committed
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
1001
1002
1003
1004
1005
1006
1007
1008
1009
1010
1011
1012
1013
1014
1015
1016
1017
1018
1019
1020
1021
1022
1023
1024
1025
1026
1027
1028
1029
1030
1031
1032
1033
1034
        {
            snappedPoint[pointI] = snappedPoints.size();
            snappedPoints.append(collapsedPoint[pointI]);
        }
    }
}


Foam::triSurface Foam::isoSurface::stitchTriPoints
(
    const bool checkDuplicates,
    const List<point>& triPoints,
    labelList& triPointReverseMap,  // unmerged to merged point
    labelList& triMap               // merged to unmerged triangle
) const
{
    label nTris = triPoints.size()/3;

    if ((triPoints.size() % 3) != 0)
    {
        FatalErrorIn("isoSurface::stitchTriPoints(..)")
            << "Problem: number of points " << triPoints.size()
            << " not a multiple of 3." << abort(FatalError);
    }

    pointField newPoints;
    mergePoints
    (
        triPoints,
        mergeDistance_,
        false,
        triPointReverseMap,
        newPoints
    );

    // Check that enough merged.
    if (debug)
    {
        pointField newNewPoints;
        labelList oldToNew;
        bool hasMerged = mergePoints
        (
            newPoints,
            mergeDistance_,
            true,
            oldToNew,
            newNewPoints
        );

        if (hasMerged)
        {
            FatalErrorIn("isoSurface::stitchTriPoints(..)")
                << "Merged points contain duplicates"
                << " when merging with distance " << mergeDistance_ << endl
                << "merged:" << newPoints.size() << " re-merged:"
                << newNewPoints.size()
                << abort(FatalError);
        }
    }


    List<labelledTri> tris;
    {
        DynamicList<labelledTri> dynTris(nTris);
        label rawPointI = 0;
        DynamicList<label> newToOldTri(nTris);

        for (label oldTriI = 0; oldTriI < nTris; oldTriI++)
        {
            labelledTri tri
            (
                triPointReverseMap[rawPointI],
                triPointReverseMap[rawPointI+1],
                triPointReverseMap[rawPointI+2],
                0
            );
            rawPointI += 3;

            if ((tri[0] != tri[1]) && (tri[0] != tri[2]) && (tri[1] != tri[2]))
            {
                newToOldTri.append(oldTriI);
                dynTris.append(tri);
            }
        }

1035
1036
        triMap.transfer(newToOldTri);
        tris.transfer(dynTris);
mattijs's avatar
mattijs committed
1037
1038
    }

mattijs's avatar
mattijs committed
1039

mattijs's avatar
mattijs committed
1040

mattijs's avatar
mattijs committed
1041
    // Determine 'flat hole' situation (see RMT paper).
mattijs's avatar
mattijs committed
1042
1043
1044
1045
1046
    // Two unconnected triangles get connected because (some of) the edges
    // separating them get collapsed. Below only checks for duplicate triangles,
    // not non-manifold edge connectivity.
    if (checkDuplicates)
    {
mattijs's avatar
mattijs committed
1047
1048
1049
1050
1051
        labelListList pointFaces;
        invertManyToMany(newPoints.size(), tris, pointFaces);

        // Filter out duplicates.
        DynamicList<label> newToOldTri(tris.size());
mattijs's avatar
mattijs committed
1052
1053
1054

        forAll(tris, triI)
        {
mattijs's avatar
mattijs committed
1055
1056
            const labelledTri& tri = tris[triI];
            const labelList& pFaces = pointFaces[tri[0]];
mattijs's avatar
mattijs committed
1057

mattijs's avatar
mattijs committed
1058
1059
1060
            // Find the maximum of any duplicates. Maximum since the tris
            // below triI
            // get overwritten so we cannot use them in a comparison.
mattijs's avatar
mattijs committed
1061
1062
            label dupTriI = -1;
            forAll(pFaces, i)
mattijs's avatar
mattijs committed
1063
            {
mattijs's avatar
mattijs committed
1064
1065
1066
                label nbrTriI = pFaces[i];

                if (nbrTriI > triI && (tris[nbrTriI] == tri))
mattijs's avatar
mattijs committed
1067
                {
mattijs's avatar
mattijs committed
1068
1069
1070
                    //Pout<< "Duplicate : " << triI << " verts:" << tri
                    //    << " to " << nbrTriI << " verts:" << tris[nbrTriI]
                    //    << endl;
mattijs's avatar
mattijs committed
1071
1072
                    dupTriI = nbrTriI;
                    break;
mattijs's avatar
mattijs committed
1073
1074
1075
                }
            }

mattijs's avatar
mattijs committed
1076
1077
            if (dupTriI == -1)
            {
mattijs's avatar
mattijs committed
1078
                // There is no (higher numbered) duplicate triangle
mattijs's avatar
mattijs committed
1079
1080
1081
1082
                label newTriI = newToOldTri.size();
                newToOldTri.append(triI);
                tris[newTriI] = tris[triI];
            }
mattijs's avatar
mattijs committed
1083
        }
mattijs's avatar
mattijs committed
1084
1085
1086

        triMap.transfer(newToOldTri);
        tris.setSize(triMap.size());
mattijs's avatar
mattijs committed
1087
1088
1089
1090
1091
1092

        if (debug)
        {
            Pout<< "isoSurface : merged from " << nTris
                << " down to " << tris.size() << " unique triangles." << endl;
        }
mattijs's avatar
mattijs committed
1093
1094
1095
1096
1097
1098
1099
1100
1101
1102
1103
1104
1105
1106
1107
1108
1109
1110
1111
1112
1113
1114
1115
1116
1117
1118
1119
1120
1121
1122
1123
1124
1125
1126
1127
1128
1129


        if (debug)
        {
            triSurface surf(tris, geometricSurfacePatchList(0), newPoints);

            forAll(surf, faceI)
            {
                const labelledTri& f = surf[faceI];
                const labelList& fFaces = surf.faceFaces()[faceI];

                forAll(fFaces, i)
                {
                    label nbrFaceI = fFaces[i];

                    if (nbrFaceI <= faceI)
                    {
                        // lower numbered faces already checked
                        continue;
                    }

                    const labelledTri& nbrF = surf[nbrFaceI];

                    if (f == nbrF)
                    {
                        FatalErrorIn("validTri(const triSurface&, const label)")
                            << "Check : "
                            << " triangle " << faceI << " vertices " << f
                            << " fc:" << f.centre(surf.points())
                            << " has the same vertices as triangle " << nbrFaceI
                            << " vertices " << nbrF
                            << " fc:" << nbrF.centre(surf.points())
                            << abort(FatalError);
                    }
                }
            }
        }
mattijs's avatar
mattijs committed
1130
1131
1132
1133
1134
1135
1136
1137
1138
1139
1140
1141
1142
1143
1144
1145
1146
1147
1148
1149
1150
1151
1152
1153
1154
1155
1156
1157
1158
1159
1160
1161
1162
1163
1164
1165
1166
1167
1168
1169
1170
1171
1172
1173
1174
1175
1176
1177
1178
1179
1180
1181
1182
1183
1184
1185
1186
1187
1188
1189
1190
1191
1192
1193
    }

    return triSurface(tris, geometricSurfacePatchList(0), newPoints, true);
}


// Does face use valid vertices?
bool Foam::isoSurface::validTri(const triSurface& surf, const label faceI)
{
    // Simple check on indices ok.

    const labelledTri& f = surf[faceI];

    if
    (
        (f[0] < 0) || (f[0] >= surf.points().size())
     || (f[1] < 0) || (f[1] >= surf.points().size())
     || (f[2] < 0) || (f[2] >= surf.points().size())
    )
    {
        WarningIn("validTri(const triSurface&, const label)")
            << "triangle " << faceI << " vertices " << f
            << " uses point indices outside point range 0.."
            << surf.points().size()-1 << endl;

        return false;
    }

    if ((f[0] == f[1]) || (f[0] == f[2]) || (f[1] == f[2]))
    {
        WarningIn("validTri(const triSurface&, const label)")
            << "triangle " << faceI
            << " uses non-unique vertices " << f
            << endl;
        return false;
    }

    // duplicate triangle check

    const labelList& fFaces = surf.faceFaces()[faceI];

    // Check if faceNeighbours use same points as this face.
    // Note: discards normal information - sides of baffle are merged.
    forAll(fFaces, i)
    {
        label nbrFaceI = fFaces[i];

        if (nbrFaceI <= faceI)
        {
            // lower numbered faces already checked
            continue;
        }

        const labelledTri& nbrF = surf[nbrFaceI];

        if
        (
            ((f[0] == nbrF[0]) || (f[0] == nbrF[1]) || (f[0] == nbrF[2]))
         && ((f[1] == nbrF[0]) || (f[1] == nbrF[1]) || (f[1] == nbrF[2]))
         && ((f[2] == nbrF[0]) || (f[2] == nbrF[1]) || (f[2] == nbrF[2]))
        )
        {
            WarningIn("validTri(const triSurface&, const label)")
                << "triangle " << faceI << " vertices " << f
mattijs's avatar
mattijs committed
1194
                << " fc:" << f.centre(surf.points())
mattijs's avatar
mattijs committed
1195
1196
                << " has the same vertices as triangle " << nbrFaceI
                << " vertices " << nbrF
mattijs's avatar
mattijs committed
1197
                << " fc:" << nbrF.centre(surf.points())
mattijs's avatar
mattijs committed
1198
1199
1200
1201
1202
1203
1204
1205
1206
1207
1208
1209
1210
1211
1212
1213
1214
1215
1216
1217
1218
1219
1220
1221
1222
1223
1224
1225
1226
1227
1228
1229
1230
1231
1232
1233
1234
1235
1236
1237
1238
1239
1240
1241
1242
                << endl;

            return false;
        }
    }
    return true;
}


void Foam::isoSurface::calcAddressing
(
    const triSurface& surf,
    List<FixedList<label, 3> >& faceEdges,
    labelList& edgeFace0,
    labelList& edgeFace1,
    Map<labelList>& edgeFacesRest
) const
{
    const pointField& points = surf.points();

    pointField edgeCentres(3*surf.size());
    label edgeI = 0;
    forAll(surf, triI)
    {
        const labelledTri& tri = surf[triI];
        edgeCentres[edgeI++] = 0.5*(points[tri[0]]+points[tri[1]]);
        edgeCentres[edgeI++] = 0.5*(points[tri[1]]+points[tri[2]]);
        edgeCentres[edgeI++] = 0.5*(points[tri[2]]+points[tri[0]]);
    }

    pointField mergedCentres;
    labelList oldToMerged;
    bool hasMerged = mergePoints
    (
        edgeCentres,
        mergeDistance_,
        false,
        oldToMerged,
        mergedCentres
    );

    if (debug)
    {
        Pout<< "isoSurface : detected "
            << mergedCentres.size()
mattijs's avatar
mattijs committed
1243
            << " geometric edges on " << surf.size() << " triangles." << endl;
mattijs's avatar
mattijs committed
1244
1245
1246
1247
1248
1249
1250
1251
1252
1253
1254
1255
1256
1257
1258
1259
1260
    }

    if (!hasMerged)
    {
        return;
    }


    // Determine faceEdges
    faceEdges.setSize(surf.size());
    edgeI = 0;
    forAll(surf, triI)
    {
        faceEdges[triI][0] = oldToMerged[edgeI++];
        faceEdges[triI][1] = oldToMerged[edgeI++];
        faceEdges[triI][2] = oldToMerged[edgeI++];
    }
mattijs's avatar
mattijs committed
1261

mattijs's avatar
mattijs committed
1262
1263
1264
1265
1266
1267
1268
1269

    // Determine edgeFaces
    edgeFace0.setSize(mergedCentres.size());
    edgeFace0 = -1;
    edgeFace1.setSize(mergedCentres.size());
    edgeFace1 = -1;
    edgeFacesRest.clear();

mattijs's avatar
mattijs committed
1270
1271
1272
1273
    // Overflow edge faces for geometric shared edges that turned
    // out to be different anyway.
    EdgeMap<labelList> extraEdgeFaces(mergedCentres.size()/100);

mattijs's avatar
mattijs committed
1274
1275
1276
1277
1278
1279
1280
    forAll(oldToMerged, oldEdgeI)
    {
        label triI = oldEdgeI / 3;
        label edgeI = oldToMerged[oldEdgeI];

        if (edgeFace0[edgeI] == -1)
        {
mattijs's avatar
mattijs committed
1281
            // First triangle for edge
mattijs's avatar
mattijs committed
1282
1283
1284
1285
            edgeFace0[edgeI] = triI;
        }
        else
        {
mattijs's avatar
mattijs committed
1286
1287
1288
1289
1290
1291
1292
1293
1294
1295
1296
1297
            //- Check that the two triangles actually topologically
            //  share an edge
            const labelledTri& prevTri = surf[edgeFace0[edgeI]];
            const labelledTri& tri = surf[triI];

            label fp = oldEdgeI % 3;

            edge e(tri[fp], tri[tri.fcIndex(fp)]);

            label prevTriIndex = -1;

            forAll(prevTri, i)
mattijs's avatar
mattijs committed
1298
            {
mattijs's avatar
mattijs committed
1299
1300
1301
1302
1303
1304
1305
1306
1307
1308
1309
1310
1311
1312
1313
1314
1315
1316
1317
1318
1319
1320
1321
1322
1323
1324
1325
                if (edge(prevTri[i], prevTri[prevTri.fcIndex(i)]) == e)
                {
                    prevTriIndex = i;
                    break;
                }
            }

            if (prevTriIndex == -1)
            {
                // Different edge. Store for later.
                EdgeMap<labelList>::iterator iter = extraEdgeFaces.find(e);

                if (iter != extraEdgeFaces.end())
                {
                    labelList& eFaces = iter();
                    label sz = eFaces.size();
                    eFaces.setSize(sz+1);
                    eFaces[sz] = triI;
                }
                else
                {
                    extraEdgeFaces.insert(e, labelList(1, triI));
                }
            }
            else if (edgeFace1[edgeI] == -1)
            {
                edgeFace1[edgeI] = triI;
mattijs's avatar
mattijs committed
1326
1327
1328
            }
            else
            {
mattijs's avatar
mattijs committed
1329
1330
1331
1332
1333
1334
1335
1336
1337
1338
1339
1340
1341
1342
1343
1344
1345
1346
1347
1348
1349
1350
1351
1352
1353
1354
1355
1356
1357
1358
1359
1360
1361
1362
1363
1364
1365
1366
1367
1368
1369
1370
1371
1372
1373
1374
1375
1376
1377
1378
1379
1380
1381
1382
1383
1384
1385
1386
1387
1388
1389
1390
1391
1392
1393
1394
1395
1396
1397
                //WarningIn("orientSurface(triSurface&)")
                //    << "Edge " << edgeI << " with centre "
                //    << mergedCentres[edgeI]
                //    << " used by more than two triangles: "
                //    << edgeFace0[edgeI] << ", "
                //    << edgeFace1[edgeI] << " and " << triI << endl;
                Map<labelList>::iterator iter = edgeFacesRest.find(edgeI);

                if (iter != edgeFacesRest.end())
                {
                    labelList& eFaces = iter();
                    label sz = eFaces.size();
                    eFaces.setSize(sz+1);
                    eFaces[sz] = triI;
                }
                else
                {
                    edgeFacesRest.insert(edgeI, labelList(1, triI));
                }
            }
        }
    }

    // Add extraEdgeFaces
    edgeI = edgeFace0.size();

    edgeFace0.setSize(edgeI + extraEdgeFaces.size());
    edgeFace1.setSize(edgeI + extraEdgeFaces.size(), -1);

    forAllConstIter(EdgeMap<labelList>, extraEdgeFaces, iter)
    {
        const labelList& eFaces = iter();

        // The current edge will become edgeI. Replace all occurrences in
        // faceEdges
        forAll(eFaces, i)
        {
            label triI = eFaces[i];
            const labelledTri& tri = surf[triI];

            FixedList<label, 3>& fEdges = faceEdges[triI];
            forAll(tri, fp)
            {
                edge e(tri[fp], tri[tri.fcIndex(fp)]);

                if (e == iter.key())
                {
                    fEdges[fp] = edgeI;
                    break;
                }
            }
        }


        // Add face to edgeFaces

        edgeFace0[edgeI] = eFaces[0];

        if (eFaces.size() >= 2)
        {
            edgeFace1[edgeI] = eFaces[1];

            if (eFaces.size() > 2)
            {
                edgeFacesRest.insert
                (
                    edgeI,
                    SubList<label>(eFaces, eFaces.size()-2, 2)
                );
mattijs's avatar
mattijs committed
1398
1399
            }
        }
mattijs's avatar
mattijs committed
1400
1401

        edgeI++;
mattijs's avatar
mattijs committed
1402
1403
1404
1405
1406
1407
1408
1409
1410
1411
1412
1413
1414
1415
1416
1417
1418
1419
1420
    }
}


void Foam::isoSurface::walkOrientation
(
    const triSurface& surf,
    const List<FixedList<label, 3> >& faceEdges,
    const labelList& edgeFace0,
    const labelList& edgeFace1,
    const label seedTriI,
    labelList& flipState
)
{
    // Do walk for consistent orientation.
    DynamicList<label> changedFaces(surf.size());

    changedFaces.append(seedTriI);

1421
    while (changedFaces.size())
mattijs's avatar
mattijs committed
1422
1423
1424
1425
1426
1427
1428
1429
1430
1431
1432
1433
    {
        DynamicList<label> newChangedFaces(changedFaces.size());

        forAll(changedFaces, i)
        {
            label triI = changedFaces[i];
            const labelledTri& tri = surf[triI];
            const FixedList<label, 3>& fEdges = faceEdges[triI];

            forAll(fEdges, fp)
            {
                label edgeI = fEdges[fp];
1434

mattijs's avatar
mattijs committed
1435
1436
1437
1438
1439
1440
1441
1442
1443
1444
1445
1446
1447
1448
1449
1450
1451
                // my points:
                label p0 = tri[fp];
                label p1 = tri[tri.fcIndex(fp)];

                label nbrI =
                (
                    edgeFace0[edgeI] != triI
                  ? edgeFace0[edgeI]
                  : edgeFace1[edgeI]
                );

                if (nbrI != -1 && flipState[nbrI] == -1)
                {
                    const labelledTri& nbrTri = surf[nbrI];

                    // nbr points
                    label nbrFp = findIndex(nbrTri, p0);
mattijs's avatar
mattijs committed
1452
1453
1454
1455
1456
1457
1458
1459
1460
1461
1462
1463
1464
1465
1466
1467
1468
1469

                    if (nbrFp == -1)
                    {
                        FatalErrorIn("isoSurface::walkOrientation(..)")
                            << "triI:" << triI
                            << " tri:" << tri
                            << " p0:" << p0
                            << " p1:" << p1
                            << " fEdges:" << fEdges
                            << " edgeI:" << edgeI
                            << " edgeFace0:" << edgeFace0[edgeI]
                            << " edgeFace1:" << edgeFace1[edgeI]
                            << " nbrI:" << nbrI
                            << " nbrTri:" << nbrTri
                            << abort(FatalError);
                    }


mattijs's avatar
mattijs committed
1470
1471
1472
1473
1474
1475
1476
1477
1478
1479
1480
1481
1482
1483
1484
1485
1486
1487
1488
                    label nbrP1 = nbrTri[nbrTri.rcIndex(nbrFp)];

                    bool sameOrientation = (p1 == nbrP1);

                    if (flipState[triI] == 0)
                    {
                        flipState[nbrI] = (sameOrientation ? 0 : 1);
                    }
                    else
                    {
                        flipState[nbrI] = (sameOrientation ? 1 : 0);
                    }
                    newChangedFaces.append(nbrI);
                }
            }
        }

        changedFaces.transfer(newChangedFaces);
    }
1489
}
mattijs's avatar
mattijs committed
1490
1491
1492
1493
1494
1495
1496
1497
1498
1499
1500
1501
1502
1503
1504
1505
1506
1507
1508
1509
1510
1511
1512
1513


void Foam::isoSurface::orientSurface
(
    triSurface& surf,
    const List<FixedList<label, 3> >& faceEdges,
    const labelList& edgeFace0,
    const labelList& edgeFace1,
    const Map<labelList>& edgeFacesRest
)
{
    // -1 : unvisited
    //  0 : leave as is
    //  1 : flip
    labelList flipState(surf.size(), -1);

    label seedTriI = 0;

    while (true)
    {
        // Find first unvisited triangle
        for
        (
            ;
1514
            seedTriI < surf.size() && flipState[seedTriI] != -1;
mattijs's avatar
mattijs committed
1515
1516
1517
1518
1519
1520
1521
1522
1523
1524
1525
1526
1527
1528
1529
1530
1531
1532
1533
1534
1535
1536
1537
1538
1539
1540
1541
1542
1543
1544
1545
1546
1547
1548
1549
1550
1551
1552
            seedTriI++
        )
        {}

        if (seedTriI == surf.size())
        {
            break;
        }

        // Note: Determine orientation of seedTriI?
        // for now assume it is ok
        flipState[seedTriI] = 0;

        walkOrientation
        (
            surf,
            faceEdges,
            edgeFace0,
            edgeFace1,
            seedTriI,
            flipState
        );
    }

    // Do actual flipping
    surf.clearOut();
    forAll(surf, triI)
    {
        if (flipState[triI] == 1)
        {
            labelledTri tri(surf[triI]);

            surf[triI][0] = tri[0];
            surf[triI][1] = tri[2];
            surf[triI][2] = tri[1];
        }
        else if (flipState[triI] == -1)
        {
mattijs's avatar
mattijs committed
1553