src/foam/interpolations/MixingPlaneInterpolation/MixingPlaneInterpolationAddressing.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | foam-extend: Open Source CFD
   4    \\    /   O peration     | Version:     3.2
   5     \\  /    A nd           | Web:         http://www.foam-extend.org
   6      \\/     M anipulation  | For copyright notice see file Copyright
   7 -------------------------------------------------------------------------------
   8 License
   9     This file is part of foam-extend.
  10
  11     foam-extend 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 3 of the License, or (at your
  14     option) any later version.
  15
  16     foam-extend is distributed in the hope that it will be useful, but
  17     WITHOUT ANY WARRANTY; without even the implied warranty of
  18     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  19     General Public License for more details.
  20
  21     You should have received a copy of the GNU General Public License
  22     along with foam-extend.  If not, see <http://www.gnu.org/licenses/>.
  23
  24 Description
  25     Mixing plane class dealing with transfer of data between two
  26     primitivePatches
  27
  28 Author
  29     Martin Beaudoin, Hydro-Quebec, 2009.  All rights reserved
  30
  31 Contributor
  32     Hrvoje Jasak, Wikki Ltd.
  33
  34 \*---------------------------------------------------------------------------*/
  35
  36 #include "MixingPlaneInterpolationTemplate.H"
  37 #include "demandDrivenData.H"
  38 #include "PrimitivePatchTemplate.H"
  39 #include "IOmanip.H"
  40 #include "transform.H"
  41 #include "RodriguesRotation.H"
  42
  43 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  44
  45 namespace Foam
  46 {
  47
  48 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  49
  50 // Compute master and slave patch addressing and weighting factors
  51 //
  52 //  Addressing:
  53 //     Given the interpolation profile, find under which profile interval falls
  54 //     any given  master and slave patch faces.
  55 //
  56 //  Weighting factors:
  57 //     Find the amount of area intersection any given patch faces has with an
  58 //     imaginary ribbon spanned from the profile interval that overlaps it.
  59 //          Weighting factor == 1.0:  The face is completely overlapped
  60 //           by its ribbon,
  61 //          Weighting factor == 0.0:  The face is not overlapped by its ribbon.
  62 //          GGI weighting factors are well suited for this.
  63 //
  64 //  For every interpolation profile intervals, we got:
  65 //        'n' addresses and weighting factors for the master patch
  66 //        'm' addresses and weighting factors for the slave patch
  67 //
  68
  69 template<class MasterPatch, class SlavePatch>
  70 void MixingPlaneInterpolation<MasterPatch, SlavePatch>::calcAddressing() const
  71 {
  72     if
  73     (
  74        masterPatchToProfileAddrPtr_
  75     || masterProfileToPatchAddrPtr_
  76     || masterPatchToProfileWeightsPtr_
  77     || masterProfileToPatchWeightsPtr_
  78     || slavePatchToProfileAddrPtr_
  79     || slaveProfileToPatchAddrPtr_
  80     || slavePatchToProfileWeightsPtr_
  81     || slaveProfileToPatchWeightsPtr_
  82     )
  83     {
  84         FatalErrorIn
  85         (
  86             "void MixingPlaneInterpolation<MasterPatch, "
  87             "SlavePatch>::calcAddressing() const"
  88         )   << "Addressing already calculated"
  89             << abort(FatalError);
  90     }
  91
  92     if (debug)
  93     {
  94         InfoIn
  95         (
  96             "void MixingPlaneInterpolation<MasterPatch, "
  97             "SlavePatch>::calcAddressing() const"
  98         )   << "Creating internal GGIs: Large values for the master GGI "
  99             << "weighting factor corrections are expected."
 100             << endl;
 101     }
 102
 103     // Construct 2 GGIs in order to evaluate the interpolation weighting
 104     // factors and the addressing
 105     // - Since the profile cannot exactly match either master and slave
 106     // patches perfectly (discretisation effects),
 107     // - it is best to use a GGI in order to compute those parameters
 108
 109     // Basic GGI with no mesh transform nor separation
 110     const tensorField noTransform(0);
 111     const vectorField noTranslation(0);
 112
 113     // NB: It is important here that the "ribbon patch" be the master patch
 114     // for each GGIs.
 115     //     We have to remember that when evaluating the GGI weighting factors,
 116     //     the slave patch faces points
 117     //     are projected onto a plane defined at every single master patch
 118     //     faces. We then compute a 2D polygon
 119     //     intersection using Sutherland-Hodgman, and so on...
 120     //     It is thus important that the "ribbon patch" be the same master
 121     //     patch for both GGI because the faces
 122     //     for both masterPatchCylCoord and slavePatchCylCoord will then be
 123     //     projected onto a common geometry.
 124     //     By doing so, we hope to minimize geometrical discretisation errors
 125     //     introduced by the possibly different
 126     //     master and slave patch mesh resolution.
 127
 128     GGIInterpolation<standAlonePatch, standAlonePatch>
 129         masterCircumAvgPatchToPatch
 130         (
 131             this->mixingPlanePatch(),
 132             this->transformedMasterPatch(),
 133             noTransform,
 134             noTransform,
 135             noTranslation,
 136             0,
 137             0,
 138             true    // Scale GGI weights
 139         );
 140
 141     GGIInterpolation<standAlonePatch, standAlonePatch>
 142         slaveCircumAvgPatchToPatch
 143         (
 144             this->mixingPlanePatch(),
 145             this->transformedShadowPatch(),
 146             noTransform,
 147             noTransform,
 148             noTranslation,
 149             0,
 150             0,
 151             true    // Scale GGI weights
 152         );
 153
 154     // Memorized the GGI addressing and weighting factors
 155     //
 156     // The master/slave weights and the master/slave addressing values
 157     // from both GGI give the information necessary to compute the
 158     // fields circumferential average on the master side, and to transfer
 159     // that averaged value properly back to the slave side.
 160     //
 161     // We get, for each master/slave patches: which face contribute to which
 162     // ribbon, and in which proportion.
 163     //
 164     // The GGI weighting factors will be use to compute the circumferential
 165     // weighted average.
 166     // Since the GGI weighting factors are already factoring in the
 167     // ratio of intersected area versus the full facet area,
 168     // we end up with a surface-area weighted average when using the GGI
 169     // weighting factors... Pretty neat... MB :)
 170
 171     // Once we collect this information, we will can simply discard the GGIs
 172     // we no longer need them
 173
 174     // The transfer from patch to profile requires the master
 175     //  weighting factors so we can do a proper facet surface weighted average
 176     masterPatchToProfileAddrPtr_ =
 177         new labelListList(masterCircumAvgPatchToPatch.masterAddr());
 178
 179     masterPatchToProfileWeightsPtr_ =
 180         new scalarListList(masterCircumAvgPatchToPatch.masterWeights());
 181
 182     slavePatchToProfileAddrPtr_ =
 183         new labelListList(slaveCircumAvgPatchToPatch.masterAddr());
 184
 185     slavePatchToProfileWeightsPtr_ =
 186         new scalarListList(slaveCircumAvgPatchToPatch.masterWeights());
 187
 188
 189     // The transfer from profile to patch requires the slave weighting
 190     // factors so we can transfer the information back to the
 191     // slave patch using the standard GGI mechanism
 192     masterProfileToPatchAddrPtr_ =
 193         new labelListList(masterCircumAvgPatchToPatch.slaveAddr());
 194
 195     masterProfileToPatchWeightsPtr_ =
 196         new scalarListList(masterCircumAvgPatchToPatch.slaveWeights());
 197
 198     slaveProfileToPatchAddrPtr_ =
 199         new labelListList(slaveCircumAvgPatchToPatch.slaveAddr());
 200
 201     slaveProfileToPatchWeightsPtr_ =
 202         new scalarListList(slaveCircumAvgPatchToPatch.slaveWeights());
 203 }
 204
 205
 206 //  Generate Cartesian/Cylindrical tranformation tensor fields for
 207 //  master and slave patches
 208 template<class MasterPatch, class SlavePatch>
 209 void
 210 MixingPlaneInterpolation<MasterPatch, SlavePatch>::calcTransforms() const
 211 {
 212     // Collapse in sweepAxis-wise direction
 213     const direction sweepDir = sweepAxisSwitch();
 214
 215     // Master side
 216     masterPatchToProfileTPtr_ = new tensorField(masterPatch_.size());
 217     tensorField& mPatchToProfileT = *masterPatchToProfileTPtr_;
 218
 219     masterProfileToPatchTPtr_ = new tensorField(masterPatch_.size());
 220     tensorField& mProfileToPatchT = *masterProfileToPatchTPtr_;
 221
 222     if (cs_.type() == coordinateSystem::typeName)
 223     {
 224         // Identity tensor for cartesian space
 225         mPatchToProfileT = tensor(sphericalTensor::I);
 226         mProfileToPatchT = tensor(sphericalTensor::I);
 227     }
 228     else
 229     {
 230         // Get master patch face centers.
 231         // Warning: Face normals are not a good choice here. We need a vector
 232         // that will intersect the rotation axis. Furthermore, if the surface
 233         // normals are all parallel with the rotation axis, no valid rotation
 234         // tensors can be computed
 235         // Use face cell centers cause the fields are taken from cell centers
 236         vectorField globalMasterVectors =
 237             masterPatch_.faceCentres() - cs_.origin();
 238
 239         // We need unit vectors for computing rotation tensor
 240         // We also need vectors lying into the plane normal to the rotation
 241         // axis. This is a major limitation of the current implementation of
 242         // rotationTensor() that is being used for computing rotation tensors:
 243         // we cannot specify the rotation axis. Instead, when using
 244         // rotationTensor(), the rotation axis will be aligned with the normal
 245         // of the plane spanned by the two vectors.
 246         // RodriguesRotation() was designed to overcome these limitations.
 247
 248         // Move face vector into the local coordinate system
 249         vectorField localMasterVectors = cs_.localVector(globalMasterVectors);
 250
 251         // Translate everything to theta=0
 252         localMasterVectors.replace(sweepDir, 0);
 253
 254         // Transform back to global
 255         vectorField transformMasterVectors =
 256             cs_.globalVector(localMasterVectors);
 257
 258         // Calculate transform tensors. These are pure rotation tensors,
 259         // aligned with the reference frame rotation axis
 260         mPatchToProfileT =
 261             RodriguesRotation
 262             (
 263                 cs_.axis(),
 264                 globalMasterVectors,
 265                 transformMasterVectors
 266             );
 267
 268         mProfileToPatchT =
 269             RodriguesRotation
 270             (
 271                 cs_.axis(),
 272                 transformMasterVectors,
 273                 globalMasterVectors
 274             );
 275     }
 276
 277     // Slave side
 278
 279     slavePatchToProfileTPtr_ = new tensorField(slavePatch_.size());
 280     tensorField& sPatchToProfileT = *slavePatchToProfileTPtr_;
 281
 282     slaveProfileToPatchTPtr_ = new tensorField(slavePatch_.size());
 283     tensorField& sProfileToPatchT = *slaveProfileToPatchTPtr_;
 284
 285     if (cs_.type() == coordinateSystem::typeName)
 286     {
 287         // Identity tensor for cartesian space
 288         sPatchToProfileT = tensor(sphericalTensor::I);
 289         sProfileToPatchT = tensor(sphericalTensor::I);
 290     }
 291     else
 292     {
 293         // Get slave patch face centers.
 294         // Warning: Face normals are not a good choice here. We need a vector
 295         // that will intersect the rotation axis. Furthermore, if the surface
 296         // normals are all parallel with the rotation axis, no valid rotation
 297         // tensors can be computed
 298         vectorField globalSlaveVectors =
 299             slavePatch_.faceCentres() - cs_.origin();
 300
 301         // We need unit vectors for computing rotation tensor
 302         // We also need vectors lying into the plane normal to the rotation
 303         // axis. This is a major limitation of the current implementation of
 304         // rotationTensor() that is being used for computing rotation tensors:
 305         // we cannot specify the rotation axis. Instead, when using
 306         // rotationTensor(), the rotation axis will be aligned with the normal
 307         // of the plane spanned by the two vectors.
 308         // RodriguesRotation() was designed to overcome these limitations.
 309
 310         // Move face vector into the local coordinate system
 311         vectorField localSlaveVectors = cs_.localVector(globalSlaveVectors);
 312
 313         // Translate everything to theta = 0
 314         localSlaveVectors.replace(sweepDir, 0);
 315
 316         // Transform back to global
 317         vectorField transformSlaveVectors =
 318             cs_.globalVector(localSlaveVectors);
 319
 320         // Calculate transform tensors. These are pure rotation tensors,
 321         // aligned with the reference frame rotation axis
 322         sPatchToProfileT =
 323             RodriguesRotation
 324             (
 325                 cs_.axis(),
 326                 globalSlaveVectors,
 327                 transformSlaveVectors
 328             );
 329
 330         sProfileToPatchT =
 331             RodriguesRotation
 332             (
 333                 cs_.axis(),
 334                 transformSlaveVectors,
 335                 globalSlaveVectors
 336             );
 337
 338         if (debug)
 339         {
 340             // We can remove this later when agree this is ok.
 341             const vectorField& slaveFaceCntr = slavePatch_.faceCentres();
 342
 343             pointField transformedFaceCentre1 = transform
 344             (
 345                 sPatchToProfileT,
 346                 slaveFaceCntr
 347             );
 348
 349             pointField transformedFaceCentre2 = transform
 350             (
 351                 sProfileToPatchT,
 352                 transformedFaceCentre1
 353             );
 354
 355             InfoIn
 356             (
 357                 "MixingPlaneInterpolation"
 358                 "<SlavePatch, SlavePatch>::calcTransforms() const"
 359             )   << "slave face centre transformation check: "
 360                 << "(should be zero!) mag = "
 361                 << mag(transformedFaceCentre2 - slaveFaceCntr) << nl
 362                 << " sum mag= "
 363                 << sum(mag(transformedFaceCentre2 - slaveFaceCntr))
 364                 << endl;
 365         }
 366     }
 367 }
 368
 369
 370 template<class MasterPatch, class SlavePatch>
 371 void
 372 MixingPlaneInterpolation<MasterPatch, SlavePatch>::clearAddressing()
 373 {
 374     // Master side
 375     deleteDemandDrivenData(masterPatchToProfileAddrPtr_);
 376     deleteDemandDrivenData(masterProfileToPatchAddrPtr_);
 377     deleteDemandDrivenData(masterPatchToProfileWeightsPtr_);
 378     deleteDemandDrivenData(masterProfileToPatchWeightsPtr_);
 379
 380     // Slave side
 381     deleteDemandDrivenData(slavePatchToProfileAddrPtr_);
 382     deleteDemandDrivenData(slaveProfileToPatchAddrPtr_);
 383     deleteDemandDrivenData(slavePatchToProfileWeightsPtr_);
 384     deleteDemandDrivenData(slaveProfileToPatchWeightsPtr_);
 385 }
 386
 387
 388 template<class MasterPatch, class SlavePatch>
 389 void
 390 MixingPlaneInterpolation<MasterPatch, SlavePatch>::clearTransforms()
 391 {
 392     // Master side
 393     deleteDemandDrivenData(masterPatchToProfileTPtr_);
 394     deleteDemandDrivenData(masterProfileToPatchTPtr_);
 395
 396     // Slave side
 397     deleteDemandDrivenData(slavePatchToProfileTPtr_);
 398     deleteDemandDrivenData(slaveProfileToPatchTPtr_);
 399 }
 400
 401
 402 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 403
 404 } // End namespace Foam
 405
 406 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 407
 408
 409 // ************************************************************************* //