polyMesh.C 32.2 KB
Newer Older
1
2
3
4
/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
5
    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
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.
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/>.
23
24
25
26
27
28
29
30
31
32

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

#include "polyMesh.H"
#include "Time.H"
#include "cellIOList.H"
#include "wedgePolyPatch.H"
#include "emptyPolyPatch.H"
#include "globalMeshData.H"
#include "processorPolyPatch.H"
33
#include "polyMeshTetDecomposition.H"
34
35
#include "indexedOctree.H"
#include "treeDataCell.H"
36
#include "MeshObject.H"
37
#include "pointMesh.H"
mattijs's avatar
mattijs committed
38

39
40
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //

41
42
namespace Foam
{
43
    defineTypeNameAndDebug(polyMesh, 0);
44

45
46
    word polyMesh::defaultRegion = "region0";
    word polyMesh::meshSubDir = "polyMesh";
47
}
48
49
50
51
52
53
54
55


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

void Foam::polyMesh::calcDirections() const
{
    for (direction cmpt=0; cmpt<vector::nComponents; cmpt++)
    {
mattijs's avatar
mattijs committed
56
        solutionD_[cmpt] = 1;
57
58
    }

mattijs's avatar
mattijs committed
59
60
61
    // Knock out empty and wedge directions. Note:they will be present on all
    // domains.

62
    label nEmptyPatches = 0;
mattijs's avatar
mattijs committed
63
    label nWedgePatches = 0;
64

mattijs's avatar
mattijs committed
65
66
    vector emptyDirVec = vector::zero;
    vector wedgeDirVec = vector::zero;
67
68
69

    forAll(boundaryMesh(), patchi)
    {
mattijs's avatar
mattijs committed
70
        if (boundaryMesh()[patchi].size())
71
        {
mattijs's avatar
mattijs committed
72
            if (isA<emptyPolyPatch>(boundaryMesh()[patchi]))
73
74
            {
                nEmptyPatches++;
mattijs's avatar
mattijs committed
75
76
77
78
79
80
81
82
83
84
85
                emptyDirVec += sum(cmptMag(boundaryMesh()[patchi].faceAreas()));
            }
            else if (isA<wedgePolyPatch>(boundaryMesh()[patchi]))
            {
                const wedgePolyPatch& wpp = refCast<const wedgePolyPatch>
                (
                    boundaryMesh()[patchi]
                );

                nWedgePatches++;
                wedgeDirVec += cmptMag(wpp.centreNormal());
86
87
88
89
            }
        }
    }

90
91
92
    reduce(nEmptyPatches, maxOp<label>());
    reduce(nWedgePatches, maxOp<label>());

93
94
    if (nEmptyPatches)
    {
mattijs's avatar
mattijs committed
95
        reduce(emptyDirVec, sumOp<vector>());
96

mattijs's avatar
mattijs committed
97
        emptyDirVec /= mag(emptyDirVec);
98
99
100

        for (direction cmpt=0; cmpt<vector::nComponents; cmpt++)
        {
mattijs's avatar
mattijs committed
101
            if (emptyDirVec[cmpt] > 1e-6)
102
            {
mattijs's avatar
mattijs committed
103
                solutionD_[cmpt] = -1;
104
105
106
            }
            else
            {
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
129
130
131
                solutionD_[cmpt] = 1;
            }
        }
    }


    // Knock out wedge directions

    geometricD_ = solutionD_;

    if (nWedgePatches)
    {
        reduce(wedgeDirVec, sumOp<vector>());

        wedgeDirVec /= mag(wedgeDirVec);

        for (direction cmpt=0; cmpt<vector::nComponents; cmpt++)
        {
            if (wedgeDirVec[cmpt] > 1e-6)
            {
                geometricD_[cmpt] = -1;
            }
            else
            {
                geometricD_[cmpt] = 1;
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
            }
        }
    }
}


// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //

Foam::polyMesh::polyMesh(const IOobject& io)
:
    objectRegistry(io),
    primitiveMesh(),
    points_
    (
        IOobject
        (
            "points",
            time().findInstance(meshDir(), "points"),
            meshSubDir,
            *this,
            IOobject::MUST_READ,
            IOobject::NO_WRITE
        )
    ),
    faces_
    (
        IOobject
        (
            "faces",
            time().findInstance(meshDir(), "faces"),
            meshSubDir,
            *this,
            IOobject::MUST_READ,
            IOobject::NO_WRITE
        )
    ),
    owner_
    (
        IOobject
        (
            "owner",
173
            faces_.instance(),
174
175
176
177
178
179
180
181
182
183
184
            meshSubDir,
            *this,
            IOobject::READ_IF_PRESENT,
            IOobject::NO_WRITE
        )
    ),
    neighbour_
    (
        IOobject
        (
            "neighbour",
185
            faces_.instance(),
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
            meshSubDir,
            *this,
            IOobject::READ_IF_PRESENT,
            IOobject::NO_WRITE
        )
    ),
    clearedPrimitives_(false),
    boundary_
    (
        IOobject
        (
            "boundary",
            time().findInstance(meshDir(), "boundary"),
            meshSubDir,
            *this,
            IOobject::MUST_READ,
            IOobject::NO_WRITE
        ),
        *this
    ),
    bounds_(points_),
207
    comm_(UPstream::worldComm),
mattijs's avatar
mattijs committed
208
209
    geometricD_(Vector<label>::zero),
    solutionD_(Vector<label>::zero),
210
    tetBasePtIsPtr_(NULL),
211
    cellTreePtr_(NULL),
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
    pointZones_
    (
        IOobject
        (
            "pointZones",
            time().findInstance
            (
                meshDir(),
                "pointZones",
                IOobject::READ_IF_PRESENT
            ),
            meshSubDir,
            *this,
            IOobject::READ_IF_PRESENT,
            IOobject::NO_WRITE
        ),
        *this
    ),
    faceZones_
    (
        IOobject
        (
            "faceZones",
            time().findInstance
            (
                meshDir(),
                "faceZones",
                IOobject::READ_IF_PRESENT
            ),
            meshSubDir,
            *this,
            IOobject::READ_IF_PRESENT,
            IOobject::NO_WRITE
        ),
        *this
    ),
    cellZones_
    (
        IOobject
        (
            "cellZones",
            time().findInstance
            (
                meshDir(),
                "cellZones",
                IOobject::READ_IF_PRESENT
            ),
            meshSubDir,
            *this,
            IOobject::READ_IF_PRESENT,
            IOobject::NO_WRITE
        ),
        *this
    ),
    globalMeshDataPtr_(NULL),
    moving_(false),
mattijs's avatar
mattijs committed
268
    topoChanging_(false),
269
270
271
272
273
274
275
276
277
    curMotionTimeIndex_(time().timeIndex()),
    oldPointsPtr_(NULL)
{
    if (exists(owner_.objectPath()))
    {
        initMesh();
    }
    else
    {
278
        cellCompactIOList cLst
279
280
281
282
283
284
285
286
287
288
289
290
291
        (
            IOobject
            (
                "cells",
                time().findInstance(meshDir(), "cells"),
                meshSubDir,
                *this,
                IOobject::MUST_READ,
                IOobject::NO_WRITE
            )
        );

        // Set the primitive mesh
292
        initMesh(cLst);
293
294
295
296
297
298
299
300
301
302
303

        owner_.write();
        neighbour_.write();
    }

    // Calculate topology for the patches (processor-processor comms etc.)
    boundary_.updateMesh();

    // Calculate the geometry for the patches (transformation tensors etc.)
    boundary_.calcGeometry();

304
305
    // Warn if global empty mesh
    if (returnReduce(nPoints(), sumOp<label>()) == 0)
306
    {
307
        WarningInFunction
308
309
            << "no points in mesh" << endl;
    }
310
    if (returnReduce(nCells(), sumOp<label>()) == 0)
311
    {
312
        WarningInFunction
313
314
            << "no cells in mesh" << endl;
    }
315
316
317

    // Initialise demand-driven data
    calcDirections();
318
319
320
321
322
323
}


Foam::polyMesh::polyMesh
(
    const IOobject& io,
Mark Olesen's avatar
Mark Olesen committed
324
325
326
327
    const Xfer<pointField>& points,
    const Xfer<faceList>& faces,
    const Xfer<labelList>& owner,
    const Xfer<labelList>& neighbour,
328
329
330
331
332
333
334
335
336
337
338
339
340
    const bool syncPar
)
:
    objectRegistry(io),
    primitiveMesh(),
    points_
    (
        IOobject
        (
            "points",
            instance(),
            meshSubDir,
            *this,
341
            io.readOpt(),
342
343
344
345
346
347
348
349
350
351
352
353
            IOobject::AUTO_WRITE
        ),
        points
    ),
    faces_
    (
        IOobject
        (
            "faces",
            instance(),
            meshSubDir,
            *this,
354
            io.readOpt(),
355
356
357
358
359
360
361
362
363
364
365
366
            IOobject::AUTO_WRITE
        ),
        faces
    ),
    owner_
    (
        IOobject
        (
            "owner",
            instance(),
            meshSubDir,
            *this,
367
            io.readOpt(),
368
369
370
371
372
373
374
375
376
377
378
379
            IOobject::AUTO_WRITE
        ),
        owner
    ),
    neighbour_
    (
        IOobject
        (
            "neighbour",
            instance(),
            meshSubDir,
            *this,
380
            io.readOpt(),
381
382
383
384
385
386
387
388
389
390
391
392
393
            IOobject::AUTO_WRITE
        ),
        neighbour
    ),
    clearedPrimitives_(false),
    boundary_
    (
        IOobject
        (
            "boundary",
            instance(),
            meshSubDir,
            *this,
394
            io.readOpt(),
395
396
397
            IOobject::AUTO_WRITE
        ),
        *this,
398
        polyPatchList()
399
400
    ),
    bounds_(points_, syncPar),
401
    comm_(UPstream::worldComm),
mattijs's avatar
mattijs committed
402
403
    geometricD_(Vector<label>::zero),
    solutionD_(Vector<label>::zero),
404
    tetBasePtIsPtr_(NULL),
405
    cellTreePtr_(NULL),
406
407
408
409
410
411
412
413
    pointZones_
    (
        IOobject
        (
            "pointZones",
            instance(),
            meshSubDir,
            *this,
414
            io.readOpt(),
415
416
417
            IOobject::NO_WRITE
        ),
        *this,
418
        PtrList<pointZone>()
419
420
421
422
423
424
425
426
427
    ),
    faceZones_
    (
        IOobject
        (
            "faceZones",
            instance(),
            meshSubDir,
            *this,
428
            io.readOpt(),
429
430
431
            IOobject::NO_WRITE
        ),
        *this,
432
        PtrList<faceZone>()
433
434
435
436
437
438
439
440
441
    ),
    cellZones_
    (
        IOobject
        (
            "cellZones",
            instance(),
            meshSubDir,
            *this,
442
            io.readOpt(),
443
444
445
            IOobject::NO_WRITE
        ),
        *this,
446
        PtrList<cellZone>()
447
448
449
    ),
    globalMeshDataPtr_(NULL),
    moving_(false),
mattijs's avatar
mattijs committed
450
    topoChanging_(false),
451
452
453
454
    curMotionTimeIndex_(time().timeIndex()),
    oldPointsPtr_(NULL)
{
    // Check if the faces and cells are valid
455
    forAll(faces_, facei)
456
    {
457
        const face& curFace = faces_[facei];
458
459
460

        if (min(curFace) < 0 || max(curFace) > points_.size())
        {
461
462
            FatalErrorInFunction
                << "Face " << facei << "contains vertex labels out of range: "
463
464
465
466
467
468
469
470
471
472
473
474
475
                << curFace << " Max point index = " << points_.size()
                << abort(FatalError);
        }
    }

    // Set the primitive mesh
    initMesh();
}


Foam::polyMesh::polyMesh
(
    const IOobject& io,
Mark Olesen's avatar
Mark Olesen committed
476
477
478
    const Xfer<pointField>& points,
    const Xfer<faceList>& faces,
    const Xfer<cellList>& cells,
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
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
    const bool syncPar
)
:
    objectRegistry(io),
    primitiveMesh(),
    points_
    (
        IOobject
        (
            "points",
            instance(),
            meshSubDir,
            *this,
            IOobject::NO_READ,
            IOobject::AUTO_WRITE
        ),
        points
    ),
    faces_
    (
        IOobject
        (
            "faces",
            instance(),
            meshSubDir,
            *this,
            IOobject::NO_READ,
            IOobject::AUTO_WRITE
        ),
        faces
    ),
    owner_
    (
        IOobject
        (
            "owner",
            instance(),
            meshSubDir,
            *this,
            IOobject::NO_READ,
            IOobject::AUTO_WRITE
        ),
        0
    ),
    neighbour_
    (
        IOobject
        (
            "neighbour",
            instance(),
            meshSubDir,
            *this,
            IOobject::NO_READ,
            IOobject::AUTO_WRITE
        ),
        0
    ),
    clearedPrimitives_(false),
    boundary_
    (
        IOobject
        (
            "boundary",
            instance(),
            meshSubDir,
            *this,
            IOobject::NO_READ,
            IOobject::AUTO_WRITE
        ),
        *this,
        0
    ),
    bounds_(points_, syncPar),
552
    comm_(UPstream::worldComm),
mattijs's avatar
mattijs committed
553
554
    geometricD_(Vector<label>::zero),
    solutionD_(Vector<label>::zero),
555
    tetBasePtIsPtr_(NULL),
556
    cellTreePtr_(NULL),
557
558
559
560
561
562
563
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
593
594
595
596
597
598
599
600
    pointZones_
    (
        IOobject
        (
            "pointZones",
            instance(),
            meshSubDir,
            *this,
            IOobject::NO_READ,
            IOobject::NO_WRITE
        ),
        *this,
        0
    ),
    faceZones_
    (
        IOobject
        (
            "faceZones",
            instance(),
            meshSubDir,
            *this,
            IOobject::NO_READ,
            IOobject::NO_WRITE
        ),
        *this,
        0
    ),
    cellZones_
    (
        IOobject
        (
            "cellZones",
            instance(),
            meshSubDir,
            *this,
            IOobject::NO_READ,
            IOobject::NO_WRITE
        ),
        *this,
        0
    ),
    globalMeshDataPtr_(NULL),
    moving_(false),
mattijs's avatar
mattijs committed
601
    topoChanging_(false),
602
603
604
    curMotionTimeIndex_(time().timeIndex()),
    oldPointsPtr_(NULL)
{
mattijs's avatar
mattijs committed
605
    // Check if faces are valid
606
    forAll(faces_, facei)
607
    {
608
        const face& curFace = faces_[facei];
609
610
611

        if (min(curFace) < 0 || max(curFace) > points_.size())
        {
612
613
            FatalErrorInFunction
                << "Face " << facei << "contains vertex labels out of range: "
614
615
616
617
618
                << curFace << " Max point index = " << points_.size()
                << abort(FatalError);
        }
    }

619
620
    // transfer in cell list
    cellList cLst(cells);
mattijs's avatar
mattijs committed
621
622

    // Check if cells are valid
623
    forAll(cLst, celli)
624
    {
625
        const cell& curCell = cLst[celli];
626
627
628

        if (min(curCell) < 0 || max(curCell) > faces_.size())
        {
629
630
            FatalErrorInFunction
                << "Cell " << celli << "contains face labels out of range: "
631
632
633
634
635
636
                << curCell << " Max face index = " << faces_.size()
                << abort(FatalError);
        }
    }

    // Set the primitive mesh
637
    initMesh(cLst);
638
639
640
641
642
}


void Foam::polyMesh::resetPrimitives
(
Mark Olesen's avatar
Mark Olesen committed
643
644
645
646
    const Xfer<pointField>& points,
    const Xfer<faceList>& faces,
    const Xfer<labelList>& owner,
    const Xfer<labelList>& neighbour,
647
648
649
650
651
    const labelList& patchSizes,
    const labelList& patchStarts,
    const bool validBoundary
)
{
mattijs's avatar
mattijs committed
652
653
    // Clear addressing. Keep geometric props and updateable props for mapping.
    clearAddressing(true);
654

mattijs's avatar
mattijs committed
655
    // Take over new primitive data.
656
    // Optimized to avoid overwriting data at all
657
    if (notNull(points))
658
    {
659
        points_.transfer(points());
660
661
        bounds_ = boundBox(points_, validBoundary);
    }
662

663
    if (notNull(faces))
664
    {
665
        faces_.transfer(faces());
666
    }
667

668
    if (notNull(owner))
669
    {
670
        owner_.transfer(owner());
671
    }
672

673
    if (notNull(neighbour))
674
    {
675
        neighbour_.transfer(neighbour());
676
677
    }

mattijs's avatar
mattijs committed
678

679
680
681
682
683
    // Reset patch sizes and starts
    forAll(boundary_, patchI)
    {
        boundary_[patchI] = polyPatch
        (
684
685
            boundary_[patchI],
            boundary_,
686
            patchI,
687
688
            patchSizes[patchI],
            patchStarts[patchI]
689
690
691
692
693
694
695
696
        );
    }


    // Flags the mesh files as being changed
    setInstance(time().timeName());

    // Check if the faces and cells are valid
697
    forAll(faces_, facei)
698
    {
699
        const face& curFace = faces_[facei];
700
701
702

        if (min(curFace) < 0 || max(curFace) > points_.size())
        {
703
704
            FatalErrorInFunction
                << "Face " << facei << " contains vertex labels out of range: "
705
706
707
708
709
710
                << curFace << " Max point index = " << points_.size()
                << abort(FatalError);
        }
    }


711
712
    // Set the primitive mesh from the owner_, neighbour_.
    // Works out from patch end where the active faces stop.
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
    initMesh();


    if (validBoundary)
    {
        // Note that we assume that all the patches stay the same and are
        // correct etc. so we can already use the patches to do
        // processor-processor comms.

        // Calculate topology for the patches (processor-processor comms etc.)
        boundary_.updateMesh();

        // Calculate the geometry for the patches (transformation tensors etc.)
        boundary_.calcGeometry();

728
729
730
731
732
733
        // Warn if global empty mesh
        if
        (
            (returnReduce(nPoints(), sumOp<label>()) == 0)
         || (returnReduce(nCells(), sumOp<label>()) == 0)
        )
734
        {
735
736
            FatalErrorInFunction
                << "no points or no cells in mesh" << endl;
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
        }
    }
}


// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //

Foam::polyMesh::~polyMesh()
{
    clearOut();
    resetMotion();
}


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

const Foam::fileName& Foam::polyMesh::dbDir() const
{
    if (objectRegistry::dbDir() == defaultRegion)
    {
        return parent().dbDir();
    }
    else
    {
        return objectRegistry::dbDir();
    }
}


Foam::fileName Foam::polyMesh::meshDir() const
{
    return dbDir()/meshSubDir;
}


const Foam::fileName& Foam::polyMesh::pointsInstance() const
{
    return points_.instance();
}


const Foam::fileName& Foam::polyMesh::facesInstance() const
{
    return faces_.instance();
}


mattijs's avatar
mattijs committed
784
const Foam::Vector<Foam::label>& Foam::polyMesh::geometricD() const
785
{
mattijs's avatar
mattijs committed
786
    if (geometricD_.x() == 0)
787
788
789
790
    {
        calcDirections();
    }

mattijs's avatar
mattijs committed
791
    return geometricD_;
792
793
794
795
796
}


Foam::label Foam::polyMesh::nGeometricD() const
{
mattijs's avatar
mattijs committed
797
798
    return cmptSum(geometricD() + Vector<label>::one)/2;
}
799
800


mattijs's avatar
mattijs committed
801
802
803
const Foam::Vector<Foam::label>& Foam::polyMesh::solutionD() const
{
    if (solutionD_.x() == 0)
804
    {
mattijs's avatar
mattijs committed
805
        calcDirections();
806
807
    }

mattijs's avatar
mattijs committed
808
    return solutionD_;
809
810
811
812
813
}


Foam::label Foam::polyMesh::nSolutionD() const
{
mattijs's avatar
mattijs committed
814
    return cmptSum(solutionD() + Vector<label>::one)/2;
815
816
817
}


818
819
const Foam::labelList& Foam::polyMesh::tetBasePtIs() const
{
820
    if (tetBasePtIsPtr_.empty())
821
822
823
    {
        if (debug)
        {
824
            WarningInFunction
825
826
827
828
                << "Forcing storage of base points."
                << endl;
        }

829
        tetBasePtIsPtr_.reset
830
        (
831
832
833
834
            new labelList
            (
                polyMeshTetDecomposition::findFaceBasePts(*this)
            )
835
836
837
        );
    }

838
839
840
841
842
843
844
845
846
847
848
849
850
    return tetBasePtIsPtr_();
}


const Foam::indexedOctree<Foam::treeDataCell>&
Foam::polyMesh::cellTree() const
{
    if (cellTreePtr_.empty())
    {
        treeBoundBox overallBb(points());

        Random rndGen(261782);

851
        overallBb = overallBb.extend(rndGen, 1e-4);
852
853
854
855
856
857
858
859
860
        overallBb.min() -= point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
        overallBb.max() += point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);

        cellTreePtr_.reset
        (
            new indexedOctree<treeDataCell>
            (
                treeDataCell
                (
861
                    false,      // not cache bb
862
                    *this,
863
                    CELL_TETS   // use tet-decomposition for any inside test
864
865
866
867
868
869
870
871
872
873
                ),
                overallBb,
                8,              // maxLevel
                10,             // leafsize
                5.0             // duplicity
            )
        );
    }

    return cellTreePtr_();
874
875
876
}


877
878
879
880
881
882
void Foam::polyMesh::addPatches
(
    const List<polyPatch*>& p,
    const bool validBoundary
)
{
883
    if (boundaryMesh().size())
884
    {
885
886
        FatalErrorInFunction
            << "boundary already exists"
887
888
889
            << abort(FatalError);
    }

mattijs's avatar
mattijs committed
890
891
892
893
    // Reset valid directions
    geometricD_ = Vector<label>::zero;
    solutionD_ = Vector<label>::zero;

894
895
896
    boundary_.setSize(p.size());

    // Copy the patch pointers
897
    forAll(p, pI)
898
899
900
901
902
903
904
905
    {
        boundary_.set(pI, p[pI]);
    }

    // parallelData depends on the processorPatch ordering so force
    // recalculation. Problem: should really be done in removeBoundary but
    // there is some info in parallelData which might be interesting inbetween
    // removeBoundary and addPatches.
906
    globalMeshDataPtr_.clear();
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927

    if (validBoundary)
    {
        // Calculate topology for the patches (processor-processor comms etc.)
        boundary_.updateMesh();

        // Calculate the geometry for the patches (transformation tensors etc.)
        boundary_.calcGeometry();

        boundary_.checkDefinition();
    }
}


void Foam::polyMesh::addZones
(
    const List<pointZone*>& pz,
    const List<faceZone*>& fz,
    const List<cellZone*>& cz
)
{
928
    if (pointZones().size() || faceZones().size() || cellZones().size())
929
    {
930
931
        FatalErrorInFunction
            << "point, face or cell zone already exists"
932
933
934
935
936
937
938
939
940
            << abort(FatalError);
    }

    // Point zones
    if (pz.size())
    {
        pointZones_.setSize(pz.size());

        // Copy the zone pointers
941
        forAll(pz, pI)
942
943
944
945
946
947
948
949
950
951
952
953
954
        {
            pointZones_.set(pI, pz[pI]);
        }

        pointZones_.writeOpt() = IOobject::AUTO_WRITE;
    }

    // Face zones
    if (fz.size())
    {
        faceZones_.setSize(fz.size());

        // Copy the zone pointers
955
        forAll(fz, fI)
956
957
958
959
960
961
962
963
964
965
966
967
968
        {
            faceZones_.set(fI, fz[fI]);
        }

        faceZones_.writeOpt() = IOobject::AUTO_WRITE;
    }

    // Cell zones
    if (cz.size())
    {
        cellZones_.setSize(cz.size());

        // Copy the zone pointers
969
        forAll(cz, cI)
970
971
972
973
974
975
976
977
978
979
980
981
982
        {
            cellZones_.set(cI, cz[cI]);
        }

        cellZones_.writeOpt() = IOobject::AUTO_WRITE;
    }
}


const Foam::pointField& Foam::polyMesh::points() const
{
    if (clearedPrimitives_)
    {
983
        FatalErrorInFunction
984
985
986
987
988
989
990
991
            << "points deallocated"
            << abort(FatalError);
    }

    return points_;
}


992
993
994
995
996
997
998
999
1000
1001
1002
1003
bool Foam::polyMesh::upToDatePoints(const regIOobject& io) const
{
    return io.upToDate(points_);
}


void Foam::polyMesh::setUpToDatePoints(regIOobject& io) const
{
    io.eventNo() = points_.eventNo();
}


1004
1005
1006
1007
const Foam::faceList& Foam::polyMesh::faces() const
{
    if (clearedPrimitives_)
    {
1008
        FatalErrorInFunction
1009
1010
1011
1012
1013
1014
1015
1016
1017
1018
1019
1020
1021
1022
1023
1024
1025
1026
1027
1028
1029
1030
            << "faces deallocated"
            << abort(FatalError);
    }

    return faces_;
}


const Foam::labelList& Foam::polyMesh::faceOwner() const
{
    return owner_;
}


const Foam::labelList& Foam::polyMesh::faceNeighbour() const
{
    return neighbour_;
}


const Foam::pointField& Foam::polyMesh::oldPoints() const
{
1031
    if (oldPointsPtr_.empty())
1032
1033
1034
    {
        if (debug)
        {
1035
            WarningInFunction
1036
1037
1038
                << endl;
        }

1039
        oldPointsPtr_.reset(new pointField(points_));
1040
1041
1042
        curMotionTimeIndex_ = time().timeIndex();
    }

1043
    return oldPointsPtr_();
1044
1045
1046
1047
1048
}


Foam::tmp<Foam::scalarField> Foam::polyMesh::movePoints
(
1049
    const pointField& newPoints
1050
1051
1052
1053
)
{
    if (debug)
    {
1054
1055
        InfoInFunction
            << "Moving points for time " << time().value()
1056
1057
1058
1059
1060
1061
1062
1063
1064
            << " index " << time().timeIndex() << endl;
    }

    moving(true);

    // Pick up old points
    if (curMotionTimeIndex_ != time().timeIndex())
    {
        // Mesh motion in the new time step
1065
1066
        oldPointsPtr_.clear();
        oldPointsPtr_.reset(new pointField(points_));
1067
1068
1069
1070
1071
        curMotionTimeIndex_ = time().timeIndex();
    }

    points_ = newPoints;

andy's avatar
andy committed
1072
    bool moveError = false;
1073
1074
1075
    if (debug)
    {
        // Check mesh motion
1076
        if (checkMeshMotion(points_, true))
1077
        {
andy's avatar
andy committed
1078
1079
            moveError = true;

1080
            InfoInFunction
1081
1082
1083
1084
1085
1086
1087
1088
                << "Moving the mesh with given points will "
                << "invalidate the mesh." << nl
                << "Mesh motion should not be executed." << endl;
        }
    }

    points_.writeOpt() = IOobject::AUTO_WRITE;
    points_.instance() = time().timeName();
1089
    points_.eventNo() = getEvent();
1090
1091
1092
1093
1094
1095
1096
1097

    tmp<scalarField> sweptVols = primitiveMesh::movePoints
    (
        points_,
        oldPoints()
    );

    // Adjust parallel shared points
1098
    if (globalMeshDataPtr_.valid())
1099
    {
1100
        globalMeshDataPtr_().movePoints(points_);
1101
1102
1103
1104
1105
1106
1107
1108
1109
1110
1111
    }

    // Force recalculation of all geometric data with new points

    bounds_ = boundBox(points_);
    boundary_.movePoints(points_);

    pointZones_.movePoints(points_);
    faceZones_.movePoints(points_);
    cellZones_.movePoints(points_);

mattijs's avatar
mattijs committed
1112
1113
1114
1115
    // Reset valid directions (could change with rotation)
    geometricD_ = Vector<label>::zero;
    solutionD_ = Vector<label>::zero;

1116
    meshObject::movePoints<polyMesh>(*this);
1117
    meshObject::movePoints<pointMesh>(*this);
mattijs's avatar
mattijs committed
1118

1119
1120
    const_cast<Time&>(time()).functionObjects().movePoints(*this);

andy's avatar
andy committed
1121
1122
1123

    if (debug && moveError)
    {
1124
1125
1126
        // Write mesh to ease debugging. Note we want to avoid calling
        // e.g. fvMesh::write since meshPhi not yet complete.
        polyMesh::write();
andy's avatar
andy committed
1127
1128
    }

1129
1130
1131
1132
1133
1134
1135
    return sweptVols;
}


void Foam::polyMesh::resetMotion() const
{
    curMotionTimeIndex_ = 0;
1136
    oldPointsPtr_.clear();
1137
1138
1139
1140
1141
}


const Foam::globalMeshData& Foam::polyMesh::globalData() const
{
1142
    if (globalMeshDataPtr_.empty())
1143
1144
1145
1146
    {
        if (debug)
        {
            Pout<< "polyMesh::globalData() const : "
mattijs's avatar
mattijs committed
1147
                << "Constructing parallelData from processor topology"
1148
1149
1150
                << endl;
        }
        // Construct globalMeshData using processorPatch information only.
1151
        globalMeshDataPtr_.reset(new globalMeshData(*this));
1152
1153
    }

1154
    return globalMeshDataPtr_();
1155
1156
1157
}


1158
1159
1160
1161
1162
1163
1164
1165
1166
1167
1168
1169
Foam::label Foam::polyMesh::comm() const
{
    return comm_;
}


Foam::label& Foam::polyMesh::comm()
{
    return comm_;
}


1170
1171
void Foam::polyMesh::removeFiles(const fileName& instanceDir) const
{
1172
    fileName meshFilesPath = thisDb().time().path()/instanceDir/meshDir();
1173
1174
1175
1176
1177
1178
1179
1180
1181
1182
1183
1184
1185
1186

    rm(meshFilesPath/"points");
    rm(meshFilesPath/"faces");
    rm(meshFilesPath/"owner");
    rm(meshFilesPath/"neighbour");
    rm(meshFilesPath/"cells");
    rm(meshFilesPath/"boundary");
    rm(meshFilesPath/"pointZones");
    rm(meshFilesPath/"faceZones");
    rm(meshFilesPath/"cellZones");
    rm(meshFilesPath/"meshModifiers");
    rm(meshFilesPath/"parallelData");

    // remove subdirectories
1187
    if (isDir(meshFilesPath/"sets"))
1188
1189
1190
1191
1192
    {
        rmDir(meshFilesPath/"sets");
    }
}

1193

1194
1195
1196
1197
1198
1199
void Foam::polyMesh::removeFiles() const
{
    removeFiles(instance());
}


1200
1201
void Foam::polyMesh::findCellFacePt
(
1202
1203
1204
1205
    const point& p,
    label& celli,
    label& tetFacei,
    label& tetPti
1206
1207
) const
{
1208
1209
1210
    celli = -1;
    tetFacei = -1;
    tetPti = -1;
1211
1212
1213
1214

    const indexedOctree<treeDataCell>& tree = cellTree();

    // Find nearest cell to the point
1215
    pointIndexHit info = tree.findNearest(p, sqr(GREAT));
1216
1217
1218
1219
1220
1221

    if (info.hit())
    {
        label nearestCellI = tree.shapes().cellLabels()[info.index()];

        // Check the nearest cell to see if the point is inside.
1222
        findTetFacePt(nearestCellI, p, tetFacei, tetPti);
1223

1224
        if (tetFacei != -1)
1225
1226
1227
        {
            // Point was in the nearest cell

1228
            celli = nearestCellI;
1229
1230
1231
1232
1233
1234
1235

            return;
        }
        else
        {
            // Check the other possible cells that the point may be in

1236
            labelList testCells = tree.findIndices(p);
1237
1238
1239
1240
1241
1242
1243
1244
1245
1246
1247
1248
1249

            forAll(testCells, pCI)
            {
                label testCellI = tree.shapes().cellLabels()[testCells[pCI]];

                if (testCellI == nearestCellI)
                {
                    // Don't retest the nearest cell

                    continue;
                }

                // Check the test cell to see if the point is inside.
1250
                findTetFacePt(testCellI, p, tetFacei, tetPti);
1251

1252
                if (tetFacei != -1)
1253
1254
1255
                {
                    // Point was in the test cell

1256
                    celli = testCellI;
1257
1258
1259
1260
1261
1262
1263
1264

                    return;
                }
            }
        }
    }
    else
    {
1265
1266
        FatalErrorInFunction
            << "Did not find nearest cell in search tree."
1267
1268
1269
1270
1271
1272
1273
            << abort(FatalError);
    }
}


void Foam::polyMesh::findTetFacePt
(
1274
1275
1276
1277
    const label celli,
    const point& p,
    label& tetFacei,
    label& tetPti
1278
1279
1280
1281
) const
{
    const polyMesh& mesh = *this;

1282
1283
1284
    tetIndices tet(polyMeshTetDecomposition::findTet(mesh, celli, p));
    tetFacei = tet.face();
    tetPti = tet.tetPt();
1285
1286
1287
}


1288
1289
1290
bool Foam::polyMesh::pointInCell
(
    const point& p,
1291
1292
    label celli,
    const cellDecomposition decompMode
1293
1294
1295
1296
) const
{
    switch (decompMode)
    {
1297
        case FACE_PLANES:
1298
        {
1299
            return primitiveMesh::pointInCell(p, celli);
1300
1301
1302
        }
        break;

1303
        case FACE_CENTRE_TRIS:
1304
        {
1305
1306
            // only test that point is on inside of plane defined by cell face
            // triangles
1307
            const cell& cFaces = cells()[celli];
1308

1309
            forAll(cFaces, cFacei)
1310
            {
1311
1312
1313
1314
                label facei = cFaces[cFacei];
                const face& f = faces_[facei];
                const point& fc = faceCentres()[facei];
                bool isOwn = (owner_[facei] == celli);
1315
1316
1317
1318
1319
1320
1321
1322
1323
1324
1325
1326
1327
1328
1329
1330
1331

                forAll(f, fp)
                {
                    label pointI;
                    label nextPointI;

                    if (isOwn)
                    {
                        pointI = f[fp];
                        nextPointI = f.nextLabel(fp);
                    }
                    else
                    {
                        pointI = f.nextLabel(fp);
                        nextPointI = f[fp];
                    }

1332
1333
1334
1335
1336
1337
                    triPointRef faceTri
                    (
                        points()[pointI],
                        points()[nextPointI],
                        fc
                    );
1338

1339
                    vector proj = p - faceTri.centre();
1340

1341
                    if ((faceTri.normal() & proj) > 0)
1342
                    {
1343
                        return false;
1344
1345
1346
                    }
                }
            }
1347
            return true;
1348
1349
1350
        }
        break;

1351
        case FACE_DIAG_TRIS:
1352
        {
1353
1354
            // only test that point is on inside of plane defined by cell face
            // triangles
1355
            const cell& cFaces = cells()[celli];
1356

1357
            forAll(cFaces, cFacei)
1358
            {
1359
1360
                label facei = cFaces[cFacei];
                const face& f = faces_[facei];
1361

1362
                for (label tetPti = 1; tetPti < f.size() - 1; tetPti++)
1363
1364
1365
1366
1367
1368
1369
                {
                    // Get tetIndices of face triangle
                    tetIndices faceTetIs
                    (
                        polyMeshTetDecomposition::triangleTetIndices
                        (
                            *this,
1370
1371
1372
                            facei,
                            celli,
                            tetPti
1373
1374
1375
1376
1377
1378
1379
1380
1381
1382
1383
1384
1385
                        )
                    );

                    triPointRef faceTri = faceTetIs.faceTri(*this);

                    vector proj = p - faceTri.centre();

                    if ((faceTri.normal() & proj) > 0)
                    {
                        return false;
                    }
                }
            }
1386

1387
            return true;
1388
1389
        }
        break;
1390
1391
1392
1393
1394
1395
1396
1397
1398
1399
1400

        case CELL_TETS:
        {
            label tetFacei;
            label tetPti;

            findTetFacePt(celli, p, tetFacei, tetPti);

            return tetFacei != -1;
        }
        break;
1401
    }
1402

1403
1404
1405
1406
1407
1408
    return false;
}


Foam::label Foam::polyMesh::findCell
(
1409
1410
    const point& p,
    const cellDecomposition decompMode
1411
1412
) const
{
1413
1414
1415
1416
1417
    if
    (
        Pstream::parRun()
     && (decompMode == FACE_DIAG_TRIS || decompMode == CELL_TETS)
    )
1418
1419
1420
1421
1422
1423
1424
1425
1426
1427
    {
        // Force construction of face-diagonal decomposition before testing
        // for zero cells.
        //
        // If parallel running a local domain might have zero cells so never
        // construct the face-diagonal decomposition which uses parallel
        // transfers.
        (void)tetBasePtIs();
    }

1428
1429
1430
1431
1432
    if (nCells() == 0)
    {
        return -1;
    }

1433
    if (decompMode == CELL_TETS)
1434
    {
1435
1436
1437
1438
1439
1440
1441
1442
1443
1444
        // Advanced search method utilizing an octree
        // and tet-decomposition of the cells

        label celli;
        label tetFacei;
        label tetPti;

        findCellFacePt(p, celli, tetFacei, tetPti);

        return celli;
1445
    }
1446
    else
1447
    {
1448
1449
1450
1451
1452
1453
1454
1455
        // Approximate search avoiding the construction of an octree
        // and cell decomposition

        // Find the nearest cell centre to this location
        label celli = findNearestCell(p);

        // If point is in the nearest cell return
        if (pointInCell(p, celli, decompMode))
1456
        {
1457
1458
1459
1460
1461
1462
1463
            return celli;
        }
        else
        {
            // Point is not in the nearest cell so search all cells

            for (label celli = 0; celli < nCells(); celli++)
1464
            {
1465
1466
1467
1468
                if (pointInCell(p, celli, decompMode))
                {
                    return celli;
                }
1469
            }
1470
1471

            return -1;
1472
1473
1474
1475
1476
        }
    }
}


1477
// ************************************************************************* //