src/dynamicMesh/motionSmoother/motionSmootherTemplates.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 1991-2009 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 "motionSmoother.H"
  28 #include "meshTools.H"
  29 #include "processorPointPatchFields.H"
  30 #include "globalPointPatchFields.H"
  31 #include "pointConstraint.H"
  32 #include "syncTools.H"
  33
  34 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  35
  36 template<class Type>
  37 void Foam::motionSmoother::checkConstraints
  38 (
  39     GeometricField<Type, pointPatchField, pointMesh>& pf
  40 )
  41 {
  42     typedef GeometricField<Type, pointPatchField, pointMesh> FldType;
  43
  44     const polyMesh& mesh = pf.mesh();
  45
  46     const polyBoundaryMesh& bm = mesh.boundaryMesh();
  47
  48     // first count the total number of patch-patch points
  49
  50     label nPatchPatchPoints = 0;
  51
  52     forAll(bm, patchi)
  53     {
  54         if(!isA<emptyPolyPatch>(bm[patchi]))
  55         {
  56             nPatchPatchPoints += bm[patchi].boundaryPoints().size();
  57         }
  58     }
  59
  60
  61     typename FldType::GeometricBoundaryField& bFld = pf.boundaryField();
  62
  63
  64     // Evaluate in reverse order
  65
  66     forAllReverse(bFld, patchi)
  67     {
  68         bFld[patchi].initEvaluate(Pstream::blocking);   // buffered
  69     }
  70
  71     forAllReverse(bFld, patchi)
  72     {
  73         bFld[patchi].evaluate(Pstream::blocking);
  74     }
  75
  76
  77     // Save the values
  78
  79     Field<Type> boundaryPointValues(nPatchPatchPoints);
  80     nPatchPatchPoints = 0;
  81
  82     forAll(bm, patchi)
  83     {
  84         if(!isA<emptyPolyPatch>(bm[patchi]))
  85         {
  86             const labelList& bp = bm[patchi].boundaryPoints();
  87             const labelList& meshPoints = bm[patchi].meshPoints();
  88
  89             forAll(bp, pointi)
  90             {
  91                 label ppp = meshPoints[bp[pointi]];
  92                 boundaryPointValues[nPatchPatchPoints++] = pf[ppp];
  93             }
  94         }
  95     }
  96
  97
  98     // Forward evaluation
  99
 100     bFld.evaluate();
 101
 102
 103     // Check
 104
 105     nPatchPatchPoints = 0;
 106
 107     forAll(bm, patchi)
 108     {
 109         if(!isA<emptyPolyPatch>(bm[patchi]))
 110         {
 111             const labelList& bp = bm[patchi].boundaryPoints();
 112             const labelList& meshPoints = bm[patchi].meshPoints();
 113
 114             forAll(bp, pointi)
 115             {
 116                 label ppp = meshPoints[bp[pointi]];
 117
 118                 const Type& savedVal = boundaryPointValues[nPatchPatchPoints++];
 119
 120                 if (savedVal != pf[ppp])
 121                 {
 122                     FatalErrorIn
 123                     (
 124                         "motionSmoother::checkConstraints"
 125                         "(GeometricField<Type, pointPatchField, pointMesh>&)"
 126                     )   << "Patch fields are not consistent on mesh point "
 127                         << ppp << " coordinate " << mesh.points()[ppp]
 128                         << " at patch " << bm[patchi].name() << '.'
 129                         << endl
 130                         << "Reverse evaluation gives value " << savedVal
 131                         << " , forward evaluation gives value " << pf[ppp]
 132                         << abort(FatalError);
 133                 }
 134             }
 135         }
 136     }
 137 }
 138
 139
 140 template<class Type>
 141 void Foam::motionSmoother::applyCornerConstraints
 142 (
 143     GeometricField<Type, pointPatchField, pointMesh>& pf
 144 ) const
 145 {
 146     forAll(patchPatchPointConstraintPoints_, pointi)
 147     {
 148         pf[patchPatchPointConstraintPoints_[pointi]] = transform
 149         (
 150             patchPatchPointConstraintTensors_[pointi],
 151             pf[patchPatchPointConstraintPoints_[pointi]]
 152         );
 153     }
 154 }
 155
 156
 157 // Average of connected points.
 158 template <class Type>
 159 Foam::tmp<Foam::GeometricField<Type, Foam::pointPatchField, Foam::pointMesh> >
 160  Foam::motionSmoother::avg
 161 (
 162     const GeometricField<Type, pointPatchField, pointMesh>& fld,
 163     const scalarField& edgeWeight,
 164     const bool separation
 165 ) const
 166 {
 167     tmp<GeometricField<Type, pointPatchField, pointMesh> > tres
 168     (
 169         new GeometricField<Type, pointPatchField, pointMesh>
 170         (
 171             IOobject
 172             (
 173                 "avg("+fld.name()+')',
 174                 fld.time().timeName(),
 175                 fld.db(),
 176                 IOobject::NO_READ,
 177                 IOobject::NO_WRITE
 178             ),
 179             fld.mesh(),
 180             dimensioned<Type>("zero", fld.dimensions(), pTraits<Type>::zero)
 181         )
 182     );
 183     GeometricField<Type, pointPatchField, pointMesh>& res = tres();
 184
 185     const polyMesh& mesh = fld.mesh()();
 186
 187
 188     // Sum local weighted values and weights
 189     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 190
 191     // Note: on coupled edges use only one edge (through isMasterEdge)
 192     // This is done so coupled edges do not get counted double.
 193
 194     scalarField sumWeight(mesh.nPoints(), 0.0);
 195
 196     const edgeList& edges = mesh.edges();
 197
 198     forAll(edges, edgeI)
 199     {
 200         if (isMasterEdge_.get(edgeI) == 1)
 201         {
 202             const edge& e = edges[edgeI];
 203             const scalar w = edgeWeight[edgeI];
 204
 205             res[e[0]] += w*fld[e[1]];
 206             sumWeight[e[0]] += w;
 207
 208             res[e[1]] += w*fld[e[0]];
 209             sumWeight[e[1]] += w;
 210         }
 211     }
 212
 213
 214     // Add coupled contributions
 215     // ~~~~~~~~~~~~~~~~~~~~~~~~~
 216
 217     syncTools::syncPointList
 218     (
 219         mesh,
 220         res,
 221         plusEqOp<Type>(),
 222         pTraits<Type>::zero,    // null value
 223         separation              // separation
 224     );
 225
 226     syncTools::syncPointList
 227     (
 228         mesh,
 229         sumWeight,
 230         plusEqOp<scalar>(),
 231         scalar(0),              // null value
 232         false                   // separation
 233     );
 234
 235
 236     // Average
 237     // ~~~~~~~
 238
 239     forAll(res, pointI)
 240     {
 241         if (mag(sumWeight[pointI]) < VSMALL)
 242         {
 243             // Unconnected point. Take over original value
 244             res[pointI] = fld[pointI];
 245         }
 246         else
 247         {
 248             res[pointI] /= sumWeight[pointI];
 249         }
 250     }
 251
 252     res.correctBoundaryConditions();
 253     applyCornerConstraints(res);
 254
 255     return tres;
 256 }
 257
 258
 259 // smooth field (point-jacobi)
 260 template <class Type>
 261 void Foam::motionSmoother::smooth
 262 (
 263     const GeometricField<Type, pointPatchField, pointMesh>& fld,
 264     const scalarField& edgeWeight,
 265     const bool separation,
 266     GeometricField<Type, pointPatchField, pointMesh>& newFld
 267 ) const
 268 {
 269     tmp<pointVectorField> tavgFld = avg
 270     (
 271         fld,
 272         edgeWeight, // weighting
 273         separation  // whether to apply separation vector
 274     );
 275     const pointVectorField& avgFld = tavgFld();
 276
 277     forAll(fld, pointI)
 278     {
 279         if (isInternalPoint(pointI))
 280         {
 281             newFld[pointI] = 0.5*fld[pointI] + 0.5*avgFld[pointI];
 282         }
 283     }
 284
 285     newFld.correctBoundaryConditions();
 286     applyCornerConstraints(newFld);
 287 }
 288
 289
 290 //- Test synchronisation of pointField
 291 template<class Type, class CombineOp>
 292 void Foam::motionSmoother::testSyncField
 293 (
 294     const Field<Type>& fld,
 295     const CombineOp& cop,
 296     const Type& zero,
 297     const bool separation
 298 ) const
 299 {
 300     if (debug)
 301     {
 302         Pout<< "testSyncField : testing synchronisation of Field<Type>."
 303             << endl;
 304     }
 305
 306     Field<Type> syncedFld(fld);
 307
 308     syncTools::syncPointList
 309     (
 310         mesh_,
 311         syncedFld,
 312         cop,
 313         zero,       // null value
 314         separation  // separation
 315     );
 316
 317     forAll(syncedFld, i)
 318     {
 319         if (syncedFld[i] != fld[i])
 320         {
 321             FatalErrorIn
 322             (
 323                 "motionSmoother::testSyncField"
 324                 "(const Field<Type>&, const CombineOp&"
 325                 ", const Type&, const bool)"
 326             )   << "On element " << i << " value:" << fld[i]
 327                 << " synchronised value:" << syncedFld[i]
 328                 << abort(FatalError);
 329         }
 330     }
 331 }
 332
 333
 334 // ************************************************************************* //