Merge branch 'master' of ssh://git.code.sf.net/p/foam-extend/foam-extend-3.2
[foam-extend-3.2.git] / src / foam / interpolations / MixingPlaneInterpolation / MixingPlaneInterpolationAddressing.C
blobc92f5150360e04e8744ad4897d11c5d0fc219d97
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.
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.
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.
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/>.
24 Description
25     Mixing plane class dealing with transfer of data between two
26     primitivePatches
28 Author
29     Martin Beaudoin, Hydro-Quebec, 2009.  All rights reserved
31 Contributor
32     Hrvoje Jasak, Wikki Ltd.
34 \*---------------------------------------------------------------------------*/
36 #include "MixingPlaneInterpolationTemplate.H"
37 #include "demandDrivenData.H"
38 #include "PrimitivePatchTemplate.H"
39 #include "IOmanip.H"
40 #include "transform.H"
41 #include "RodriguesRotation.H"
43 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
45 namespace Foam
48 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
50 // Compute master and slave patch addressing and weighting factors
52 //  Addressing:
53 //     Given the interpolation profile, find under which profile interval falls
54 //     any given  master and slave patch faces.
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.
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
69 template<class MasterPatch, class SlavePatch>
70 void MixingPlaneInterpolation<MasterPatch, SlavePatch>::calcAddressing() const
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     }
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     }
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
109     // Basic GGI with no mesh transform nor separation
110     const tensorField noTransform(0);
111     const vectorField noTranslation(0);
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.
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         );
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         );
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 :)
171     // Once we collect this information, we will can simply discard the GGIs
172     // we no longer need them
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());
179     masterPatchToProfileWeightsPtr_ =
180         new scalarListList(masterCircumAvgPatchToPatch.masterWeights());
182     slavePatchToProfileAddrPtr_ =
183         new labelListList(slaveCircumAvgPatchToPatch.masterAddr());
185     slavePatchToProfileWeightsPtr_ =
186         new scalarListList(slaveCircumAvgPatchToPatch.masterWeights());
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());
195     masterProfileToPatchWeightsPtr_ =
196         new scalarListList(masterCircumAvgPatchToPatch.slaveWeights());
198     slaveProfileToPatchAddrPtr_ =
199         new labelListList(slaveCircumAvgPatchToPatch.slaveAddr());
201     slaveProfileToPatchWeightsPtr_ =
202         new scalarListList(slaveCircumAvgPatchToPatch.slaveWeights());
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
212     // Collapse in sweepAxis-wise direction
213     const direction sweepDir = sweepAxisSwitch();
215     // Master side
216     masterPatchToProfileTPtr_ = new tensorField(masterPatch_.size());
217     tensorField& mPatchToProfileT = *masterPatchToProfileTPtr_;
219     masterProfileToPatchTPtr_ = new tensorField(masterPatch_.size());
220     tensorField& mProfileToPatchT = *masterProfileToPatchTPtr_;
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();
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.
248         // Move face vector into the local coordinate system
249         vectorField localMasterVectors = cs_.localVector(globalMasterVectors);
251         // Translate everything to theta=0
252         localMasterVectors.replace(sweepDir, 0);
254         // Transform back to global
255         vectorField transformMasterVectors =
256             cs_.globalVector(localMasterVectors);
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             );
268         mProfileToPatchT =
269             RodriguesRotation
270             (
271                 cs_.axis(),
272                 transformMasterVectors,
273                 globalMasterVectors
274             );
275     }
277     // Slave side
279     slavePatchToProfileTPtr_ = new tensorField(slavePatch_.size());
280     tensorField& sPatchToProfileT = *slavePatchToProfileTPtr_;
282     slaveProfileToPatchTPtr_ = new tensorField(slavePatch_.size());
283     tensorField& sProfileToPatchT = *slaveProfileToPatchTPtr_;
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();
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.
310         // Move face vector into the local coordinate system
311         vectorField localSlaveVectors = cs_.localVector(globalSlaveVectors);
313         // Translate everything to theta = 0
314         localSlaveVectors.replace(sweepDir, 0);
316         // Transform back to global
317         vectorField transformSlaveVectors =
318             cs_.globalVector(localSlaveVectors);
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             );
330         sProfileToPatchT =
331             RodriguesRotation
332             (
333                 cs_.axis(),
334                 transformSlaveVectors,
335                 globalSlaveVectors
336             );
338         if (debug)
339         {
340             // We can remove this later when agree this is ok.
341             const vectorField& slaveFaceCntr = slavePatch_.faceCentres();
343             pointField transformedFaceCentre1 = transform
344             (
345                 sPatchToProfileT,
346                 slaveFaceCntr
347             );
349             pointField transformedFaceCentre2 = transform
350             (
351                 sProfileToPatchT,
352                 transformedFaceCentre1
353             );
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     }
370 template<class MasterPatch, class SlavePatch>
371 void
372 MixingPlaneInterpolation<MasterPatch, SlavePatch>::clearAddressing()
374     // Master side
375     deleteDemandDrivenData(masterPatchToProfileAddrPtr_);
376     deleteDemandDrivenData(masterProfileToPatchAddrPtr_);
377     deleteDemandDrivenData(masterPatchToProfileWeightsPtr_);
378     deleteDemandDrivenData(masterProfileToPatchWeightsPtr_);
380     // Slave side
381     deleteDemandDrivenData(slavePatchToProfileAddrPtr_);
382     deleteDemandDrivenData(slaveProfileToPatchAddrPtr_);
383     deleteDemandDrivenData(slavePatchToProfileWeightsPtr_);
384     deleteDemandDrivenData(slaveProfileToPatchWeightsPtr_);
388 template<class MasterPatch, class SlavePatch>
389 void
390 MixingPlaneInterpolation<MasterPatch, SlavePatch>::clearTransforms()
392     // Master side
393     deleteDemandDrivenData(masterPatchToProfileTPtr_);
394     deleteDemandDrivenData(masterProfileToPatchTPtr_);
396     // Slave side
397     deleteDemandDrivenData(slavePatchToProfileTPtr_);
398     deleteDemandDrivenData(slaveProfileToPatchTPtr_);
402 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
404 } // End namespace Foam
406 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
409 // ************************************************************************* //