src/dynamicMesh/dynamicFvMesh/dynamicTopoFvMesh/lengthScaleEstimator/lengthScaleEstimatorI.H

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