src/dynamicMesh/dynamicTopoFvMesh/lengthScaleEstimator/lengthScaleEstimatorI.H

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | foam-extend: Open Source CFD
   4    \\    /   O peration     | Version:     3.2
   5     \\  /    A nd           | Web:         http://www.foam-extend.org
   6      \\/     M anipulation  | For copyright notice see file Copyright
   7 -------------------------------------------------------------------------------
   8 License
   9     This file is part of foam-extend.
  10
  11     foam-extend is free software: you can redistribute it and/or modify it
  12     under the terms of the GNU General Public License as published by the
  13     Free Software Foundation, either version 3 of the License, or (at your
  14     option) any later version.
  15
  16     foam-extend is distributed in the hope that it will be useful, but
  17     WITHOUT ANY WARRANTY; without even the implied warranty of
  18     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  19     General Public License for more details.
  20
  21     You should have received a copy of the GNU General Public License
  22     along with foam-extend.  If not, see <http://www.gnu.org/licenses/>.
  23
  24 \*---------------------------------------------------------------------------*/
  25
  26 namespace Foam
  27 {
  28
  29 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  30
  31
  32 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
  33
  34 inline scalar lengthScaleEstimator::ratioMin() const
  35 {
  36     return ratioMin_;
  37 }
  38
  39
  40 inline scalar lengthScaleEstimator::ratioMax() const
  41 {
  42     return ratioMax_;
  43 }
  44
  45
  46 inline scalar lengthScaleEstimator::growthFactor() const
  47 {
  48     return growthFactor_;
  49 }
  50
  51
  52 //- Limit length scale for surface-edges
  53 inline void lengthScaleEstimator::limitScale(scalar& scale) const
  54 {
  55     // Need to account for multiplication by ratioMin / ratioMax as well
  56     scale = Foam::max(scale, minLengthScale_ / ratioMin_);
  57     scale = Foam::min(scale, maxLengthScale_ / ratioMax_);
  58 }
  59
  60
  61 //- Check if a particular patch is free-floating
  62 //  (i.e., gets its length-scale from the interior)
  63 inline bool lengthScaleEstimator::isFreePatch
  64 (
  65     const label pIndex
  66 ) const
  67 {
  68     const polyBoundaryMesh& boundary = mesh_.boundaryMesh();
  69
  70     if (freePatches_.found(boundary[pIndex].name()))
  71     {
  72         return true;
  73     }
  74
  75     // Certain boundary types are also considered floating
  76     const word& pType = boundary[pIndex].type();
  77
  78     if
  79     (
  80         (pType == "processor") ||
  81         (pType == "cyclic") ||
  82         (pType == "wedge") ||
  83         (pType == "empty") ||
  84         (pType == "symmetryPlane")
  85     )
  86     {
  87         return true;
  88     }
  89
  90     return false;
  91 }
  92
  93
  94 //- Check if a particular patch is flagged
  95 //  for proximity-based refinement
  96 inline bool lengthScaleEstimator::isProximityPatch
  97 (
  98     const label pIndex
  99 ) const
 100 {
 101     if (proximityPatches_.found(mesh_.boundaryMesh()[pIndex].name()))
 102     {
 103         return true;
 104     }
 105
 106     return false;
 107 }
 108
 109
 110 //- Check if a particular patch is flagged
 111 //  for curvature-based refinement
 112 inline bool lengthScaleEstimator::isCurvaturePatch
 113 (
 114     const label pIndex
 115 ) const
 116 {
 117     if (curvaturePatches_.found(mesh_.boundaryMesh()[pIndex].name()))
 118     {
 119         return true;
 120     }
 121
 122     return false;
 123 }
 124
 125
 126 //- Check whether a particular point is too close
 127 //  to a previous mesh slice location
 128 inline bool lengthScaleEstimator::checkOldSlices
 129 (
 130     const vector& gCentre
 131 ) const
 132 {
 133     forAll(sliceBoxes_, boxI)
 134     {
 135         if (sliceBoxes_[boxI].contains(gCentre))
 136         {
 137             // Too close to another slice-point. Bail out.
 138             return true;
 139         }
 140     }
 141
 142     return false;
 143 }
 144
 145
 146 //- Add a boundBox to the existing set of sliceBoxes
 147 inline void lengthScaleEstimator::appendBox
 148 (
 149     const boundBox& bBox
 150 )
 151 {
 152     label currentSize = sliceBoxes_.size();
 153     sliceBoxes_.setSize(currentSize + 1);
 154     sliceBoxes_[currentSize] = bBox;
 155 }
 156
 157
 158 //- Clear the list of sliceBoxes
 159 inline void lengthScaleEstimator::clearBoxes()
 160 {
 161     sliceBoxes_.clear();
 162 }
 163
 164
 165 //- Check whether a particular patch permits refinement
 166 //  - Return true if not permitted
 167 inline bool lengthScaleEstimator::checkRefinementPatch
 168 (
 169     const label pIndex
 170 ) const
 171 {
 172     if (pIndex < 0)
 173     {
 174         return false;
 175     }
 176
 177     if (findIndex(noModPatchIDs_, pIndex) > -1)
 178     {
 179         return true;
 180     }
 181
 182     // Refinement is allowed
 183     return false;
 184 }
 185
 186
 187 // Return the fixed length-scale value for a boundary face
 188 inline scalar lengthScaleEstimator::fixedLengthScale
 189 (
 190     const label fIndex,
 191     const label pIndex,
 192     bool usePolyMesh
 193 ) const
 194 {
 195     scalar scale = 0.0;
 196
 197     const polyBoundaryMesh& boundary = mesh_.boundaryMesh();
 198
 199     // Check fixed length-scale patches
 200     // If the value is negative, average face length-scales.
 201     if (fixedPatches_.found(boundary[pIndex].name()))
 202     {
 203         scalar dictValue =
 204         (
 205             fixedPatches_[boundary[pIndex].name()][0].scalarToken()
 206         );
 207
 208         if (dictValue > 0.0)
 209         {
 210             return dictValue;
 211         }
 212     }
 213
 214     // Approximate a length-scale from face area
 215     if (usePolyMesh)
 216     {
 217         scale = Foam::sqrt(2.0 * mag(mesh_.faceAreas()[fIndex]));
 218     }
 219
 220     return scale;
 221 }
 222
 223
 224 // Test an edge / face for proximity with other faces on proximity patches
 225 // and return the scalar distance to an oppositely-oriented face.
 226 inline bool lengthScaleEstimator::testProximity
 227 (
 228     const vector& gCentre,
 229     const vector& gNormal,
 230     const scalar testStep,
 231     label& proxFace,
 232     scalar& proxDistance
 233 ) const
 234 {
 235     // Reset input
 236     proxDistance = GREAT;
 237
 238     dynamicLabelList posIndices(20);
 239     scalar minDeviation = -0.9;
 240     label nD = spatialRes_, binSize = proximityBins_.size();
 241
 242     // Obtain min/max extents
 243     const point& bMin = proxBoundBox_.min();
 244     const point& bMax = proxBoundBox_.max();
 245
 246     // Extend bounding-box dimensions a bit to avoid edge-effects.
 247     scalar ext = 0.02*(mag(bMax - bMin));
 248
 249     // Define an inverse grid-cell size.
 250     scalar xL = nD/(bMax.x() - bMin.x() + ext);
 251     scalar yL = nD/(bMax.y() - bMin.y() + ext);
 252     scalar zL = nD/(bMax.z() - bMin.z() + ext);
 253
 254     // Reset the proximity face
 255     proxFace = -1;
 256
 257     // Now take multiple steps in both normal directions,
 258     // and add to the list of boxes to be checked.
 259     for (scalar dir = -1.0; dir < 2.0; dir += 2.0)
 260     {
 261         for (scalar step = 0.0; step < 5.0*testStep; step += testStep)
 262         {
 263             // Hash the point-location
 264             point p = (gCentre + (dir*step*gNormal)) - bMin;
 265
 266             label i = label(mag(::floor(p.x()*xL)));
 267             label j = label(mag(::floor(p.y()*yL)));
 268             label k = label(mag(::floor(p.z()*zL)));
 269
 270             label pos = mag(((k*nD*nD)+(j*nD)+i) % binSize);
 271
 272             if (findIndex(posIndices, pos) == -1)
 273             {
 274                 posIndices.append(pos);
 275             }
 276         }
 277     }
 278
 279     // Obtain old-mesh face geometry for reference.
 280     const vectorField& faceAreas = mesh_.faceAreas();
 281     const vectorField& faceCentres = mesh_.faceCentres();
 282
 283     forAll(posIndices, indexI)
 284     {
 285         const labelList& posBin = proximityBins_[posIndices[indexI]];
 286
 287         forAll(posBin, faceI)
 288         {
 289             // Step 1: Measure the distance to the face.
 290             vector rFace = (faceCentres[posBin[faceI]] - gCentre);
 291
 292             scalar distance = mag(rFace);
 293
 294             // Step 2: Check if this face is oriented away from face / edge.
 295             const vector& fNormal = faceAreas[posBin[faceI]];
 296
 297             scalar deviation = (gNormal & (fNormal/(mag(fNormal) + VSMALL)));
 298
 299             if
 300             (
 301                 (deviation < minDeviation) &&
 302                 (distance < proxDistance)
 303             )
 304             {
 305                 // Update statistics
 306                 proxFace = posBin[faceI];
 307                 proxDistance = distance;
 308                 // minDeviation = deviation;
 309             }
 310         }
 311     }
 312
 313     // Check for fall below threshold
 314     if ((proxFace != -1) && (proxDistance < sliceThreshold_))
 315     {
 316         return true;
 317     }
 318
 319     // Does not fall below threshold
 320     return false;
 321 }
 322
 323
 324 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 325
 326 } // End namespace Foam
 327
 328 // ************************************************************************* //