src/fvMotionSolver/pointPatchFields/derived/surfaceSlipDisplacement/surfaceSlipDisplacementPointPatchVectorField.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 1991-2008 OpenCFD Ltd.
   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 #include "surfaceSlipDisplacementPointPatchVectorField.H"
  28 #include "addToRunTimeSelectionTable.H"
  29 #include "Time.H"
  30 #include "transformField.H"
  31 #include "fvMesh.H"
  32 #include "displacementLaplacianFvMotionSolver.H"
  33
  34 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  35
  36 namespace Foam
  37 {
  38
  39 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  40
  41 template<>
  42 const char*
  43 NamedEnum<surfaceSlipDisplacementPointPatchVectorField::followMode, 3>::
  44 names[] =
  45 {
  46     "nearest",
  47     "pointNormal",
  48     "fixedNormal"
  49 };
  50
  51 const NamedEnum<surfaceSlipDisplacementPointPatchVectorField::followMode, 3>
  52     surfaceSlipDisplacementPointPatchVectorField::followModeNames_;
  53
  54
  55 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
  56
  57 surfaceSlipDisplacementPointPatchVectorField::
  58 surfaceSlipDisplacementPointPatchVectorField
  59 (
  60     const pointPatch& p,
  61     const DimensionedField<vector, pointMesh>& iF
  62 )
  63 :
  64     pointPatchVectorField(p, iF),
  65     projectMode_(NEAREST),
  66     projectDir_(vector::zero),
  67     wedgePlane_(-1)
  68 {}
  69
  70
  71 surfaceSlipDisplacementPointPatchVectorField::
  72 surfaceSlipDisplacementPointPatchVectorField
  73 (
  74     const pointPatch& p,
  75     const DimensionedField<vector, pointMesh>& iF,
  76     const dictionary& dict
  77 )
  78 :
  79     pointPatchVectorField(p, iF, dict),
  80     surfacesDict_(dict.subDict("geometry")),
  81     projectMode_(followModeNames_.read(dict.lookup("followMode"))),
  82     projectDir_(dict.lookup("projectDirection")),
  83     wedgePlane_(readLabel(dict.lookup("wedgePlane"))),
  84     frozenPointsZone_(dict.lookup("frozenPointsZone"))
  85 {}
  86
  87
  88 surfaceSlipDisplacementPointPatchVectorField::
  89 surfaceSlipDisplacementPointPatchVectorField
  90 (
  91     const surfaceSlipDisplacementPointPatchVectorField& ppf,
  92     const pointPatch& p,
  93     const DimensionedField<vector, pointMesh>& iF,
  94     const pointPatchFieldMapper&
  95 )
  96 :
  97     pointPatchVectorField(p, iF),
  98     surfacesDict_(ppf.surfacesDict()),
  99     projectMode_(ppf.projectMode()),
 100     projectDir_(ppf.projectDir()),
 101     wedgePlane_(ppf.wedgePlane()),
 102     frozenPointsZone_(ppf.frozenPointsZone())
 103 {}
 104
 105
 106 surfaceSlipDisplacementPointPatchVectorField::
 107 surfaceSlipDisplacementPointPatchVectorField
 108 (
 109     const surfaceSlipDisplacementPointPatchVectorField& ppf
 110 )
 111 :
 112     pointPatchVectorField(ppf),
 113     surfacesDict_(ppf.surfacesDict()),
 114     projectMode_(ppf.projectMode()),
 115     projectDir_(ppf.projectDir()),
 116     wedgePlane_(ppf.wedgePlane()),
 117     frozenPointsZone_(ppf.frozenPointsZone())
 118 {}
 119
 120
 121 surfaceSlipDisplacementPointPatchVectorField::
 122 surfaceSlipDisplacementPointPatchVectorField
 123 (
 124     const surfaceSlipDisplacementPointPatchVectorField& ppf,
 125     const DimensionedField<vector, pointMesh>& iF
 126 )
 127 :
 128     pointPatchVectorField(ppf, iF),
 129     surfacesDict_(ppf.surfacesDict()),
 130     projectMode_(ppf.projectMode()),
 131     projectDir_(ppf.projectDir()),
 132     wedgePlane_(ppf.wedgePlane()),
 133     frozenPointsZone_(ppf.frozenPointsZone())
 134 {}
 135
 136
 137 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 138
 139 const searchableSurfaces& surfaceSlipDisplacementPointPatchVectorField::
 140 surfaces() const
 141 {
 142     if (!surfacesPtr_.valid())
 143     {
 144         surfacesPtr_.reset
 145         (
 146             new searchableSurfaces
 147             (
 148                 IOobject
 149                 (
 150                     "abc",                              // dummy name
 151                     db().time().constant(),             // directory
 152                     "triSurface",                       // instance
 153                     db().time(),                        // registry
 154                     IOobject::MUST_READ,
 155                     IOobject::NO_WRITE
 156                 ),
 157                 surfacesDict_
 158             )
 159         );
 160     }
 161     return surfacesPtr_();
 162 }
 163
 164
 165 void surfaceSlipDisplacementPointPatchVectorField::evaluate
 166 (
 167     const Pstream::commsTypes commsType
 168 )
 169 {
 170     const polyMesh& mesh = patch().boundaryMesh().mesh()();
 171
 172     //const scalar deltaT = mesh.time().deltaT().value();
 173
 174     // Construct large enough vector in direction of projectDir so
 175     // we're guaranteed to hit something.
 176
 177     const scalar projectLen = mag(mesh.bounds().max()-mesh.bounds().min());
 178
 179     // For case of fixed projection vector:
 180     vector projectVec;
 181     if (projectMode_ == FIXEDNORMAL)
 182     {
 183         vector n = projectDir_/mag(projectDir_);
 184         projectVec = projectLen*n;
 185     }
 186
 187     //- Per point projection vector:
 188
 189     const pointField& localPoints = patch().localPoints();
 190     const labelList& meshPoints = patch().meshPoints();
 191
 192     vectorField displacement(this->patchInternalField());
 193
 194
 195     // Get fixed points (bit of a hack)
 196     const pointZone* zonePtr = NULL;
 197
 198     if (frozenPointsZone_.size() > 0)
 199     {
 200         const pointZoneMesh& pZones = mesh.pointZones();
 201
 202         zonePtr = &pZones[pZones.findZoneID(frozenPointsZone_)];
 203
 204         Pout<< "surfaceSlipDisplacementPointPatchVectorField : Fixing all "
 205             << zonePtr->size() << " points in pointZone " << zonePtr->name()
 206             << endl;
 207     }
 208
 209     // Get the starting locations from the motionSolver
 210     const displacementLaplacianFvMotionSolver& motionSolver =
 211         mesh.lookupObject<displacementLaplacianFvMotionSolver>
 212         (
 213             "dynamicMeshDict"
 214         );
 215     const pointField& points0 = motionSolver.points0();
 216
 217
 218 //XXXXXX
 219
 220
 221     pointField start(meshPoints.size());
 222     forAll(start, i)
 223     {
 224         start[i] = points0[meshPoints[i]] + displacement[i];
 225     }
 226
 227     if (projectMode_ == NEAREST)
 228     {
 229         List<pointIndexHit> nearest;
 230         labelList hitSurfaces;
 231         surfaces().findNearest
 232         (
 233             start,
 234             scalarField(start.size(), sqr(projectLen)),
 235             hitSurfaces,
 236             nearest
 237         );
 238
 239         forAll(nearest, i)
 240         {
 241             if (zonePtr && (zonePtr->whichPoint(meshPoints[i]) >= 0))
 242             {
 243                 // Fixed point. Reset to point0 location.
 244                 displacement[i] = points0[meshPoints[i]] - localPoints[i];
 245             }
 246             else if (nearest[i].hit())
 247             {
 248                 displacement[i] =
 249                     nearest[i].hitPoint()
 250                   - points0[meshPoints[i]];
 251             }
 252             else
 253             {
 254                 Pout<< "    point:" << meshPoints[i]
 255                     << " coord:" << localPoints[i]
 256                     << "  did not find any surface within " << projectLen
 257                     << endl;
 258             }
 259         }
 260     }
 261     else
 262     {
 263         // Do tests on all points. Combine later on.
 264
 265         // 1. Check if already on surface
 266         List<pointIndexHit> nearest;
 267         {
 268             labelList nearestSurface;
 269             surfaces().findNearest
 270             (
 271                 start,
 272                 scalarField(start.size(), sqr(SMALL)),
 273                 nearestSurface,
 274                 nearest
 275             );
 276         }
 277
 278         // 2. intersection. (combined later on with information from nearest
 279         // above)
 280         vectorField projectVecs(start.size(), projectVec);
 281
 282         if (projectMode_ == POINTNORMAL)
 283         {
 284             projectVecs = projectLen*patch().pointNormals();
 285         }
 286
 287         // Knock out any wedge component
 288         scalarField offset(start.size(), 0.0);
 289         if (wedgePlane_ >= 0 && wedgePlane_ <= vector::nComponents)
 290         {
 291             forAll(offset, i)
 292             {
 293                 offset[i] = start[i][wedgePlane_];
 294                 start[i][wedgePlane_] = 0;
 295                 projectVecs[i][wedgePlane_] = 0;
 296             }
 297         }
 298
 299         List<pointIndexHit> rightHit;
 300         {
 301             labelList rightSurf;
 302             surfaces().findAnyIntersection
 303             (
 304                 start,
 305                 start+projectVecs,
 306                 rightSurf,
 307                 rightHit
 308             );
 309         }
 310
 311         List<pointIndexHit> leftHit;
 312         {
 313             labelList leftSurf;
 314             surfaces().findAnyIntersection
 315             (
 316                 start,
 317                 start-projectVecs,
 318                 leftSurf,
 319                 leftHit
 320             );
 321         }
 322
 323         // 3. Choose either -fixed, nearest, right, left.
 324         forAll(displacement, i)
 325         {
 326             if (zonePtr && (zonePtr->whichPoint(meshPoints[i]) >= 0))
 327             {
 328                 // Fixed point. Reset to point0 location.
 329                 displacement[i] = points0[meshPoints[i]] - localPoints[i];
 330             }
 331             else if (nearest[i].hit())
 332             {
 333                 // Found nearest.
 334                 displacement[i] =
 335                     nearest[i].hitPoint()
 336                   - points0[meshPoints[i]];
 337             }
 338             else
 339             {
 340                 pointIndexHit interPt;
 341
 342                 if (rightHit[i].hit())
 343                 {
 344                     if (leftHit[i].hit())
 345                     {
 346                         if
 347                         (
 348                             magSqr(rightHit[i].hitPoint()-start[i])
 349                           < magSqr(leftHit[i].hitPoint()-start[i])
 350                         )
 351                         {
 352                             interPt = rightHit[i];
 353                         }
 354                         else
 355                         {
 356                             interPt = leftHit[i];
 357                         }
 358                     }
 359                     else
 360                     {
 361                         interPt = rightHit[i];
 362                     }
 363                 }
 364                 else
 365                 {
 366                     if (leftHit[i].hit())
 367                     {
 368                         interPt = leftHit[i];
 369                     }
 370                 }
 371
 372
 373                 if (interPt.hit())
 374                 {
 375                     if (wedgePlane_ >= 0 && wedgePlane_ <= vector::nComponents)
 376                     {
 377                         interPt.rawPoint()[wedgePlane_] += offset[i];
 378                     }
 379                     displacement[i] = interPt.rawPoint()-points0[meshPoints[i]];
 380                 }
 381                 else
 382                 {
 383                     Pout<< "    point:" << meshPoints[i]
 384                         << " coord:" << localPoints[i]
 385                         << "  did not find any intersection between ray from "
 386                         << start[i]-projectVecs[i]
 387                         << " to " << start[i]+projectVecs[i]
 388                         << endl;
 389                 }
 390             }
 391         }
 392     }
 393
 394     // Get internal field to insert values into
 395     Field<vector>& iF = const_cast<Field<vector>&>(this->internalField());
 396
 397     //setInInternalField(iF, motionU);
 398     setInInternalField(iF, displacement);
 399
 400     pointPatchVectorField::evaluate(commsType);
 401 }
 402
 403
 404 void surfaceSlipDisplacementPointPatchVectorField::write(Ostream& os) const
 405 {
 406     pointPatchVectorField::write(os);
 407     os.writeKeyword("geometry") << surfacesDict_
 408         << token::END_STATEMENT << nl;
 409     os.writeKeyword("followMode") << followModeNames_[projectMode_]
 410         << token::END_STATEMENT << nl;
 411     os.writeKeyword("projectDirection") << projectDir_
 412         << token::END_STATEMENT << nl;
 413     os.writeKeyword("wedgePlane") << wedgePlane_
 414         << token::END_STATEMENT << nl;
 415     os.writeKeyword("frozenPointsZone") << frozenPointsZone_
 416         << token::END_STATEMENT << nl;
 417 }
 418
 419
 420 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 421
 422 makePointPatchTypeField
 423 (
 424     pointPatchVectorField,
 425     surfaceSlipDisplacementPointPatchVectorField
 426 );
 427
 428
 429 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 430
 431 } // End namespace Foam
 432
 433 // ************************************************************************* //