src/fvMotionSolver/pointPatchFields/derived/surfaceDisplacement/surfaceDisplacementPointPatchVectorField.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 "surfaceDisplacementPointPatchVectorField.H"
  27 #include "addToRunTimeSelectionTable.H"
  28 #include "Time.H"
  29 #include "transformField.H"
  30 #include "fvMesh.H"
  31 #include "displacementLaplacianFvMotionSolver.H"
  32
  33 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  34
  35 namespace Foam
  36 {
  37
  38
  39 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  40
  41 template<>
  42 const char*
  43 NamedEnum<surfaceDisplacementPointPatchVectorField::projectMode, 3>::
  44 names[] =
  45 {
  46     "nearest",
  47     "pointNormal",
  48     "fixedNormal"
  49 };
  50
  51 const NamedEnum<surfaceDisplacementPointPatchVectorField::projectMode, 3>
  52     surfaceDisplacementPointPatchVectorField::projectModeNames_;
  53
  54
  55 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  56
  57 void surfaceDisplacementPointPatchVectorField::calcProjection
  58 (
  59     vectorField& displacement
  60 ) const
  61 {
  62     const polyMesh& mesh = patch().boundaryMesh().mesh()();
  63     const pointField& localPoints = patch().localPoints();
  64     const labelList& meshPoints = patch().meshPoints();
  65
  66     //const scalar deltaT = mesh.time().deltaTValue();
  67
  68     // Construct large enough vector in direction of projectDir so
  69     // we're guaranteed to hit something.
  70
  71     //- Per point projection vector:
  72     const scalar projectLen = mag(mesh.bounds().max()-mesh.bounds().min());
  73
  74     // For case of fixed projection vector:
  75     vector projectVec;
  76     if (projectMode_ == FIXEDNORMAL)
  77     {
  78         vector n = projectDir_/mag(projectDir_);
  79         projectVec = projectLen*n;
  80     }
  81
  82
  83     // Get fixed points (bit of a hack)
  84     const pointZone* zonePtr = NULL;
  85
  86     if (frozenPointsZone_.size() > 0)
  87     {
  88         const pointZoneMesh& pZones = mesh.pointZones();
  89
  90         zonePtr = &pZones[frozenPointsZone_];
  91
  92         Pout<< "surfaceDisplacementPointPatchVectorField : Fixing all "
  93             << zonePtr->size() << " points in pointZone " << zonePtr->name()
  94             << endl;
  95     }
  96
  97     // Get the starting locations from the motionSolver
  98     const pointField& points0 = mesh.lookupObject<displacementFvMotionSolver>
  99     (
 100         "dynamicMeshDict"
 101     ).points0();
 102
 103
 104     pointField start(meshPoints.size());
 105     forAll(start, i)
 106     {
 107         start[i] = points0[meshPoints[i]] + displacement[i];
 108     }
 109
 110     label nNotProjected = 0;
 111
 112     if (projectMode_ == NEAREST)
 113     {
 114         List<pointIndexHit> nearest;
 115         labelList hitSurfaces;
 116         surfaces().findNearest
 117         (
 118             start,
 119             scalarField(start.size(), sqr(projectLen)),
 120             hitSurfaces,
 121             nearest
 122         );
 123
 124         forAll(nearest, i)
 125         {
 126             if (zonePtr && (zonePtr->whichPoint(meshPoints[i]) >= 0))
 127             {
 128                 // Fixed point. Reset to point0 location.
 129                 displacement[i] = points0[meshPoints[i]] - localPoints[i];
 130             }
 131             else if (nearest[i].hit())
 132             {
 133                 displacement[i] =
 134                     nearest[i].hitPoint()
 135                   - points0[meshPoints[i]];
 136             }
 137             else
 138             {
 139                 nNotProjected++;
 140
 141                 if (debug)
 142                 {
 143                     Pout<< "    point:" << meshPoints[i]
 144                         << " coord:" << localPoints[i]
 145                         << "  did not find any surface within " << projectLen
 146                         << endl;
 147                 }
 148             }
 149         }
 150     }
 151     else
 152     {
 153         // Do tests on all points. Combine later on.
 154
 155         // 1. Check if already on surface
 156         List<pointIndexHit> nearest;
 157         {
 158             labelList nearestSurface;
 159             surfaces().findNearest
 160             (
 161                 start,
 162                 scalarField(start.size(), sqr(SMALL)),
 163                 nearestSurface,
 164                 nearest
 165             );
 166         }
 167
 168         // 2. intersection. (combined later on with information from nearest
 169         // above)
 170         vectorField projectVecs(start.size(), projectVec);
 171
 172         if (projectMode_ == POINTNORMAL)
 173         {
 174             projectVecs = projectLen*patch().pointNormals();
 175         }
 176
 177         // Knock out any wedge component
 178         scalarField offset(start.size(), 0.0);
 179         if (wedgePlane_ >= 0 && wedgePlane_ <= vector::nComponents)
 180         {
 181             forAll(offset, i)
 182             {
 183                 offset[i] = start[i][wedgePlane_];
 184                 start[i][wedgePlane_] = 0;
 185                 projectVecs[i][wedgePlane_] = 0;
 186             }
 187         }
 188
 189         List<pointIndexHit> rightHit;
 190         {
 191             labelList rightSurf;
 192             surfaces().findAnyIntersection
 193             (
 194                 start,
 195                 start+projectVecs,
 196                 rightSurf,
 197                 rightHit
 198             );
 199         }
 200
 201         List<pointIndexHit> leftHit;
 202         {
 203             labelList leftSurf;
 204             surfaces().findAnyIntersection
 205             (
 206                 start,
 207                 start-projectVecs,
 208                 leftSurf,
 209                 leftHit
 210             );
 211         }
 212
 213         // 3. Choose either -fixed, nearest, right, left.
 214         forAll(displacement, i)
 215         {
 216             if (zonePtr && (zonePtr->whichPoint(meshPoints[i]) >= 0))
 217             {
 218                 // Fixed point. Reset to point0 location.
 219                 displacement[i] = points0[meshPoints[i]] - localPoints[i];
 220             }
 221             else if (nearest[i].hit())
 222             {
 223                 // Found nearest.
 224                 displacement[i] =
 225                     nearest[i].hitPoint()
 226                   - points0[meshPoints[i]];
 227             }
 228             else
 229             {
 230                 pointIndexHit interPt;
 231
 232                 if (rightHit[i].hit())
 233                 {
 234                     if (leftHit[i].hit())
 235                     {
 236                         if
 237                         (
 238                             magSqr(rightHit[i].hitPoint()-start[i])
 239                           < magSqr(leftHit[i].hitPoint()-start[i])
 240                         )
 241                         {
 242                             interPt = rightHit[i];
 243                         }
 244                         else
 245                         {
 246                             interPt = leftHit[i];
 247                         }
 248                     }
 249                     else
 250                     {
 251                         interPt = rightHit[i];
 252                     }
 253                 }
 254                 else
 255                 {
 256                     if (leftHit[i].hit())
 257                     {
 258                         interPt = leftHit[i];
 259                     }
 260                 }
 261
 262
 263                 if (interPt.hit())
 264                 {
 265                     if (wedgePlane_ >= 0 && wedgePlane_ <= vector::nComponents)
 266                     {
 267                         interPt.rawPoint()[wedgePlane_] += offset[i];
 268                     }
 269                     displacement[i] = interPt.rawPoint()-points0[meshPoints[i]];
 270                 }
 271                 else
 272                 {
 273                     nNotProjected++;
 274
 275                     if (debug)
 276                     {
 277                         Pout<< "    point:" << meshPoints[i]
 278                             << " coord:" << localPoints[i]
 279                             << "  did not find any intersection between"
 280                             << " ray from " << start[i]-projectVecs[i]
 281                             << " to " << start[i]+projectVecs[i] << endl;
 282                     }
 283                 }
 284             }
 285         }
 286     }
 287
 288     reduce(nNotProjected, sumOp<label>());
 289
 290     if (nNotProjected > 0)
 291     {
 292         Info<< "surfaceDisplacement :"
 293             << " on patch " << patch().name()
 294             << " did not project " << nNotProjected
 295             << " out of " << returnReduce(localPoints.size(), sumOp<label>())
 296             << " points." << endl;
 297     }
 298 }
 299
 300
 301 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 302
 303 surfaceDisplacementPointPatchVectorField::
 304 surfaceDisplacementPointPatchVectorField
 305 (
 306     const pointPatch& p,
 307     const DimensionedField<vector, pointMesh>& iF
 308 )
 309 :
 310     fixedValuePointPatchVectorField(p, iF),
 311     velocity_(vector::zero),
 312     projectMode_(NEAREST),
 313     projectDir_(vector::zero),
 314     wedgePlane_(-1)
 315 {}
 316
 317
 318 surfaceDisplacementPointPatchVectorField::
 319 surfaceDisplacementPointPatchVectorField
 320 (
 321     const pointPatch& p,
 322     const DimensionedField<vector, pointMesh>& iF,
 323     const dictionary& dict
 324 )
 325 :
 326     fixedValuePointPatchVectorField(p, iF, dict),
 327     velocity_(dict.lookup("velocity")),
 328     surfacesDict_(dict.subDict("geometry")),
 329     projectMode_(projectModeNames_.read(dict.lookup("projectMode"))),
 330     projectDir_(dict.lookup("projectDirection")),
 331     wedgePlane_(dict.lookupOrDefault("wedgePlane", -1)),
 332     frozenPointsZone_(dict.lookupOrDefault("frozenPointsZone", word::null))
 333 {
 334     if (velocity_.x() < 0 || velocity_.y() < 0 || velocity_.z() < 0)
 335     {
 336         FatalErrorIn
 337         (
 338             "surfaceDisplacementPointPatchVectorField::\n"
 339             "surfaceDisplacementPointPatchVectorField\n"
 340             "(\n"
 341             "    const pointPatch& p,\n"
 342             "    const DimensionedField<vector, pointMesh>& iF,\n"
 343             "    const dictionary& dict\n"
 344             ")\n"
 345         )   << "All components of velocity have to be positive : "
 346             << velocity_ << nl
 347             << "Set velocity components to a great value if no clipping"
 348             << " necessary." << exit(FatalError);
 349     }
 350 }
 351
 352
 353 surfaceDisplacementPointPatchVectorField::
 354 surfaceDisplacementPointPatchVectorField
 355 (
 356     const surfaceDisplacementPointPatchVectorField& ppf,
 357     const pointPatch& p,
 358     const DimensionedField<vector, pointMesh>& iF,
 359     const pointPatchFieldMapper& mapper
 360 )
 361 :
 362     fixedValuePointPatchVectorField(ppf, p, iF, mapper),
 363     velocity_(ppf.velocity_),
 364     surfacesDict_(ppf.surfacesDict_),
 365     projectMode_(ppf.projectMode_),
 366     projectDir_(ppf.projectDir_),
 367     wedgePlane_(ppf.wedgePlane_),
 368     frozenPointsZone_(ppf.frozenPointsZone_)
 369 {}
 370
 371
 372 surfaceDisplacementPointPatchVectorField::
 373 surfaceDisplacementPointPatchVectorField
 374 (
 375     const surfaceDisplacementPointPatchVectorField& ppf
 376 )
 377 :
 378     fixedValuePointPatchVectorField(ppf),
 379     velocity_(ppf.velocity_),
 380     surfacesDict_(ppf.surfacesDict_),
 381     projectMode_(ppf.projectMode_),
 382     projectDir_(ppf.projectDir_),
 383     wedgePlane_(ppf.wedgePlane_),
 384     frozenPointsZone_(ppf.frozenPointsZone_)
 385 {}
 386
 387
 388 surfaceDisplacementPointPatchVectorField::
 389 surfaceDisplacementPointPatchVectorField
 390 (
 391     const surfaceDisplacementPointPatchVectorField& ppf,
 392     const DimensionedField<vector, pointMesh>& iF
 393 )
 394 :
 395     fixedValuePointPatchVectorField(ppf, iF),
 396     velocity_(ppf.velocity_),
 397     surfacesDict_(ppf.surfacesDict_),
 398     projectMode_(ppf.projectMode_),
 399     projectDir_(ppf.projectDir_),
 400     wedgePlane_(ppf.wedgePlane_),
 401     frozenPointsZone_(ppf.frozenPointsZone_)
 402 {}
 403
 404
 405 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 406
 407 const searchableSurfaces&
 408 surfaceDisplacementPointPatchVectorField::surfaces() const
 409 {
 410     if (surfacesPtr_.empty())
 411     {
 412         surfacesPtr_.reset
 413         (
 414             new searchableSurfaces
 415             (
 416                 IOobject
 417                 (
 418                     "abc",                              // dummy name
 419                     db().time().constant(),             // directory
 420                     "triSurface",                       // instance
 421                     db().time(),                        // registry
 422                     IOobject::MUST_READ,
 423                     IOobject::NO_WRITE
 424                 ),
 425                 surfacesDict_
 426             )
 427         );
 428     }
 429     return surfacesPtr_();
 430 }
 431
 432
 433 void surfaceDisplacementPointPatchVectorField::updateCoeffs()
 434 {
 435     if (this->updated())
 436     {
 437         return;
 438     }
 439
 440     const polyMesh& mesh = patch().boundaryMesh().mesh()();
 441
 442     vectorField currentDisplacement(this->patchInternalField());
 443
 444     // Calculate intersections with surface w.r.t points0.
 445     vectorField displacement(currentDisplacement);
 446     calcProjection(displacement);
 447
 448     // offset wrt current displacement
 449     vectorField offset(displacement-currentDisplacement);
 450
 451     // Clip offset to maximum displacement possible: velocity*timestep
 452
 453     const scalar deltaT = mesh.time().deltaTValue();
 454     const vector clipVelocity = velocity_*deltaT;
 455
 456     forAll(displacement, i)
 457     {
 458         vector& d = offset[i];
 459
 460         for (direction cmpt = 0; cmpt < vector::nComponents; cmpt++)
 461         {
 462             if (d[cmpt] < 0)
 463             {
 464                 d[cmpt] = max(d[cmpt], -clipVelocity[cmpt]);
 465             }
 466             else
 467             {
 468                 d[cmpt] = min(d[cmpt], clipVelocity[cmpt]);
 469             }
 470         }
 471     }
 472
 473     this->operator==(currentDisplacement+offset);
 474
 475     fixedValuePointPatchVectorField::updateCoeffs();
 476 }
 477
 478
 479 void surfaceDisplacementPointPatchVectorField::write(Ostream& os) const
 480 {
 481     fixedValuePointPatchVectorField::write(os);
 482     os.writeKeyword("velocity") << velocity_
 483         << token::END_STATEMENT << nl;
 484     os.writeKeyword("geometry") << surfacesDict_
 485         << token::END_STATEMENT << nl;
 486     os.writeKeyword("projectMode") << projectModeNames_[projectMode_]
 487         << token::END_STATEMENT << nl;
 488     os.writeKeyword("projectDirection") << projectDir_
 489         << token::END_STATEMENT << nl;
 490     os.writeKeyword("wedgePlane") << wedgePlane_
 491         << token::END_STATEMENT << nl;
 492     if (frozenPointsZone_ != word::null)
 493     {
 494         os.writeKeyword("frozenPointsZone") << frozenPointsZone_
 495             << token::END_STATEMENT << nl;
 496     }
 497 }
 498
 499
 500 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 501
 502 makePointPatchTypeField
 503 (
 504     fixedValuePointPatchVectorField,
 505     surfaceDisplacementPointPatchVectorField
 506 );
 507
 508
 509 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 510
 511 } // End namespace Foam
 512
 513 // ************************************************************************* //