src/dynamicMesh/motionSmoother/motionSmootherTemplates.C

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