src/OpenFOAM/meshes/primitiveMesh/primitiveMesh.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
   6      \\/     M anipulation  |
   7 -------------------------------------------------------------------------------
   8 License
   9     This file is part of OpenFOAM.
  10
  11     OpenFOAM is free software: you can redistribute it and/or modify it
  12     under the terms of the GNU General Public License as published by
  13     the Free Software Foundation, either version 3 of the License, or
  14     (at your option) any later version.
  15
  16     OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
  17     ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  18     FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
  19     for more details.
  20
  21     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 #include "primitiveMesh.H"
  27 #include "demandDrivenData.H"
  28 #include "indexedOctree.H"
  29 #include "treeDataCell.H"
  30
  31 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  32
  33 defineTypeNameAndDebug(Foam::primitiveMesh, 0);
  34
  35
  36 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
  37
  38 Foam::primitiveMesh::primitiveMesh()
  39 :
  40     nInternalPoints_(0),    // note: points are considered ordered on empty mesh
  41     nPoints_(0),
  42     nInternal0Edges_(-1),
  43     nInternal1Edges_(-1),
  44     nInternalEdges_(-1),
  45     nEdges_(-1),
  46     nInternalFaces_(0),
  47     nFaces_(0),
  48     nCells_(0),
  49
  50     cellShapesPtr_(NULL),
  51     edgesPtr_(NULL),
  52     ccPtr_(NULL),
  53     ecPtr_(NULL),
  54     pcPtr_(NULL),
  55
  56     cfPtr_(NULL),
  57     efPtr_(NULL),
  58     pfPtr_(NULL),
  59
  60     cePtr_(NULL),
  61     fePtr_(NULL),
  62     pePtr_(NULL),
  63     ppPtr_(NULL),
  64     cpPtr_(NULL),
  65
  66     cellTreePtr_(NULL),
  67
  68     labels_(0),
  69
  70     cellCentresPtr_(NULL),
  71     faceCentresPtr_(NULL),
  72     cellVolumesPtr_(NULL),
  73     faceAreasPtr_(NULL)
  74 {}
  75
  76
  77 // Construct from components
  78 // WARNING: ASSUMES CORRECT ORDERING OF DATA.
  79 Foam::primitiveMesh::primitiveMesh
  80 (
  81     const label nPoints,
  82     const label nInternalFaces,
  83     const label nFaces,
  84     const label nCells
  85 )
  86 :
  87     nInternalPoints_(-1),
  88     nPoints_(nPoints),
  89     nEdges_(-1),
  90     nInternalFaces_(nInternalFaces),
  91     nFaces_(nFaces),
  92     nCells_(nCells),
  93
  94     cellShapesPtr_(NULL),
  95     edgesPtr_(NULL),
  96     ccPtr_(NULL),
  97     ecPtr_(NULL),
  98     pcPtr_(NULL),
  99
 100     cfPtr_(NULL),
 101     efPtr_(NULL),
 102     pfPtr_(NULL),
 103
 104     cePtr_(NULL),
 105     fePtr_(NULL),
 106     pePtr_(NULL),
 107     ppPtr_(NULL),
 108     cpPtr_(NULL),
 109
 110     cellTreePtr_(NULL),
 111
 112     labels_(0),
 113
 114     cellCentresPtr_(NULL),
 115     faceCentresPtr_(NULL),
 116     cellVolumesPtr_(NULL),
 117     faceAreasPtr_(NULL)
 118 {}
 119
 120
 121 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 122
 123 Foam::primitiveMesh::~primitiveMesh()
 124 {
 125     clearOut();
 126 }
 127
 128
 129 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 130
 131 bool Foam::primitiveMesh::calcPointOrder
 132 (
 133     label& nInternalPoints,
 134     labelList& oldToNew,
 135     const faceList& faces,
 136     const label nInternalFaces,
 137     const label nPoints
 138 )
 139 {
 140     // Internal points are points that are not used by a boundary face.
 141
 142     // Map from old to new position
 143     oldToNew.setSize(nPoints);
 144     oldToNew = -1;
 145
 146
 147     // 1. Create compact addressing for boundary points. Start off by indexing
 148     // from 0 inside oldToNew. (shifted up later on)
 149
 150     label nBoundaryPoints = 0;
 151     for (label faceI = nInternalFaces; faceI < faces.size(); faceI++)
 152     {
 153         const face& f = faces[faceI];
 154
 155         forAll(f, fp)
 156         {
 157             label pointI = f[fp];
 158
 159             if (oldToNew[pointI] == -1)
 160             {
 161                 oldToNew[pointI] = nBoundaryPoints++;
 162             }
 163         }
 164     }
 165
 166     // Now we know the number of boundary and internal points
 167
 168     nInternalPoints = nPoints - nBoundaryPoints;
 169
 170     // Move the boundary addressing up
 171     forAll(oldToNew, pointI)
 172     {
 173         if (oldToNew[pointI] != -1)
 174         {
 175             oldToNew[pointI] += nInternalPoints;
 176         }
 177     }
 178
 179
 180     // 2. Compact the internal points. Detect whether internal and boundary
 181     // points are mixed.
 182
 183     label internalPointI = 0;
 184
 185     bool ordered = true;
 186
 187     for (label faceI = 0; faceI < nInternalFaces; faceI++)
 188     {
 189         const face& f = faces[faceI];
 190
 191         forAll(f, fp)
 192         {
 193             label pointI = f[fp];
 194
 195             if (oldToNew[pointI] == -1)
 196             {
 197                 if (pointI >= nInternalPoints)
 198                 {
 199                     ordered = false;
 200                 }
 201                 oldToNew[pointI] = internalPointI++;
 202             }
 203         }
 204     }
 205
 206     return ordered;
 207 }
 208
 209
 210 void Foam::primitiveMesh::reset
 211 (
 212     const label nPoints,
 213     const label nInternalFaces,
 214     const label nFaces,
 215     const label nCells
 216 )
 217 {
 218     clearOut();
 219
 220     nPoints_ = nPoints;
 221     nEdges_ = -1;
 222     nInternal0Edges_ = -1;
 223     nInternal1Edges_ = -1;
 224     nInternalEdges_ = -1;
 225
 226     nInternalFaces_ = nInternalFaces;
 227     nFaces_ = nFaces;
 228     nCells_ = nCells;
 229
 230     // Check if points are ordered
 231     label nInternalPoints;
 232     labelList pointMap;
 233
 234     bool isOrdered = calcPointOrder
 235     (
 236         nInternalPoints,
 237         pointMap,
 238         faces(),
 239         nInternalFaces_,
 240         nPoints_
 241     );
 242
 243     if (isOrdered)
 244     {
 245         nInternalPoints_ = nInternalPoints;
 246     }
 247     else
 248     {
 249         nInternalPoints_ = -1;
 250     }
 251
 252     if (debug)
 253     {
 254         Pout<< "primitiveMesh::reset : mesh reset to"
 255             << " nInternalPoints:" << nInternalPoints_
 256             << " nPoints:" << nPoints_
 257             << " nEdges:" << nEdges_
 258             << " nInternalFaces:" << nInternalFaces_
 259             << " nFaces:" << nFaces_
 260             << " nCells:" << nCells_
 261             << endl;
 262     }
 263 }
 264
 265
 266 void Foam::primitiveMesh::reset
 267 (
 268     const label nPoints,
 269     const label nInternalFaces,
 270     const label nFaces,
 271     const label nCells,
 272     cellList& clst
 273 )
 274 {
 275     reset
 276     (
 277         nPoints,
 278         nInternalFaces,
 279         nFaces,
 280         nCells
 281     );
 282
 283     cfPtr_ = new cellList(clst, true);
 284 }
 285
 286
 287 void Foam::primitiveMesh::reset
 288 (
 289     const label nPoints,
 290     const label nInternalFaces,
 291     const label nFaces,
 292     const label nCells,
 293     const Xfer<cellList>& clst
 294 )
 295 {
 296     reset
 297     (
 298         nPoints,
 299         nInternalFaces,
 300         nFaces,
 301         nCells
 302     );
 303
 304     cfPtr_ = new cellList(clst);
 305 }
 306
 307
 308 Foam::tmp<Foam::scalarField> Foam::primitiveMesh::movePoints
 309 (
 310     const pointField& newPoints,
 311     const pointField& oldPoints
 312 )
 313 {
 314     if (newPoints.size() <  nPoints() || oldPoints.size() < nPoints())
 315     {
 316         FatalErrorIn
 317         (
 318             "primitiveMesh::movePoints(const pointField& newPoints, "
 319             "const pointField& oldPoints)"
 320         )   << "Cannot move points: size of given point list smaller "
 321             << "than the number of active points"
 322             << abort(FatalError);
 323     }
 324
 325     // Create swept volumes
 326     const faceList& f = faces();
 327
 328     tmp<scalarField> tsweptVols(new scalarField(f.size()));
 329     scalarField& sweptVols = tsweptVols();
 330
 331     forAll(f, faceI)
 332     {
 333         sweptVols[faceI] = f[faceI].sweptVol(oldPoints, newPoints);
 334     }
 335
 336     // Force recalculation of all geometric data with new points
 337     clearGeom();
 338
 339     return tsweptVols;
 340 }
 341
 342
 343 const Foam::cellShapeList& Foam::primitiveMesh::cellShapes() const
 344 {
 345     if (!cellShapesPtr_)
 346     {
 347         calcCellShapes();
 348     }
 349
 350     return *cellShapesPtr_;
 351 }
 352
 353
 354 const Foam::indexedOctree<Foam::treeDataCell>&
 355 Foam::primitiveMesh::cellTree() const
 356 {
 357     if (!cellTreePtr_)
 358     {
 359         treeBoundBox overallBb(points());
 360
 361         Random rndGen(261782);
 362
 363         overallBb = overallBb.extend(rndGen, 1E-4);
 364         overallBb.min() -= point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
 365         overallBb.max() += point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
 366
 367         cellTreePtr_ =
 368             new indexedOctree<treeDataCell>
 369             (
 370                 treeDataCell
 371                 (
 372                     false,      // not cache bb
 373                     *this
 374                 ),
 375                 overallBb,
 376                 8,              // maxLevel
 377                 10,             // leafsize
 378                 5.0             // duplicity
 379             );
 380     }
 381
 382     return *cellTreePtr_;
 383 }
 384
 385
 386 // ************************************************************************* //