fixed writing out entries in advective bc
[OpenFOAM-1.6-ext.git] / src / OpenFOAM / coordinateSystems / coordinateRotation / coordinateRotation.C
blobe30de1c0c0e4f620cc204c60f3aacee94e291636
1 /*---------------------------------------------------------------------------*\
2   =========                 |
3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
4    \\    /   O peration     |
5     \\  /    A nd           | Copyright held by original author
6      \\/     M anipulation  |
7 -------------------------------------------------------------------------------
8 License
9     This file is part of OpenFOAM.
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.
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.
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
25 \*---------------------------------------------------------------------------*/
27 #include "coordinateRotation.H"
28 #include "dictionary.H"
29 #include "addToRunTimeSelectionTable.H"
31 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
33 namespace Foam
35     defineTypeNameAndDebug(coordinateRotation, 0);
36     defineRunTimeSelectionTable(coordinateRotation, dictionary);
39 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
41 void Foam::coordinateRotation::calcTransform
43     const vector& axis1,
44     const vector& axis2,
45     const axisOrder& order
48     vector a = axis1 / mag(axis1);
49     vector b = axis2;
51     // Absorb minor non-orthogonality into axis2
52     b = b - (b & a)*a;
54     if (mag(b) < SMALL)
55     {
56         FatalErrorIn("coordinateRotation::calcTransform()")
57             << "axis1, axis2 appear co-linear: "
58             << axis1 << ", " << axis2 << endl
59             << abort(FatalError);
60     }
62     b = b / mag(b);
63     vector c = a ^ b;
65     // the global -> local transformation
66     tensor Rtr;
67     switch (order)
68     {
69         case e1e2:
70             Rtr = tensor(a, b, c);
71             break;
73         case e2e3:
74             Rtr = tensor(c, a, b);
75             break;
77         case e3e1:
78             Rtr = tensor(b, c, a);
79             break;
81         default:
82             FatalErrorIn("coordinateRotation::calcTransform()")
83                 << "programmer error" << endl
84                 << abort(FatalError);
85             // To satisfy compiler warnings
86             Rtr = tensor::zero;
87             break;
88     }
90     // the local -> global transformation
91     tensor::operator=( Rtr.T() );
95 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
97 Foam::coordinateRotation::coordinateRotation()
99     tensor(sphericalTensor::I)
103 Foam::coordinateRotation::coordinateRotation
105     const vector& axis,
106     const vector& dir
109     tensor(sphericalTensor::I)
111     calcTransform(axis, dir, e3e1);
115 Foam::coordinateRotation::coordinateRotation
117     const dictionary& dict
120     tensor(sphericalTensor::I)
122     operator=(dict);
126 Foam::coordinateRotation::coordinateRotation
128     const vector& v,
129     const scalar angle
132     tensor(sphericalTensor::I)
134     scalar c = cos(angle);
135     scalar m = mag(v);
136     scalar sm = sin(angle)/m;
138     tensor Rtr =
139         (1 - c)*v*v/m/m
140       + tensor
141         (
142             c,         -v.z()*sm, v.y()*sm,
143             v.z()*sm,  c,         -v.x()*sm,
144             -v.y()*sm, v.x()*sm,  c
145         );
147     // the local -> global transformation
148     tensor::operator=(Rtr.T());
152 // * * * * * * * * * * * * * * * * Selectors * * * * * * * * * * * * * * * * //
154 Foam::autoPtr<Foam::coordinateRotation> Foam::coordinateRotation::New
156     const dictionary& dict
159     if (debug)
160     {
161         Pout<< "coordinateRotation::New(const dictionary&) : "
162             << "constructing coordinateRotation"
163             << endl;
164     }
166     // default type is self (alias: "axes")
167     word rotType(typeName_());
168     dict.readIfPresent("type", rotType);
170     // can (must) construct base class directly
171     if (rotType == typeName_() || rotType == "axes")
172     {
173         return autoPtr<coordinateRotation>(new coordinateRotation(dict));
174     }
177     dictionaryConstructorTable::iterator cstrIter =
178         dictionaryConstructorTablePtr_->find(rotType);
180     if (cstrIter == dictionaryConstructorTablePtr_->end())
181     {
182         FatalIOErrorIn
183         (
184             "coordinateRotation::New(const dictionary&)",
185             dict
186         )   << "Unknown coordinateRotation type "
187             << rotType << nl << nl
188             << "Valid coordinateRotation types are :" <<  nl
189             << "[default: axes " << typeName_() << "]"
190             << dictionaryConstructorTablePtr_->toc()
191             << exit(FatalIOError);
192     }
194     return autoPtr<coordinateRotation>(cstrIter()(dict));
198 // * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
200 void Foam::coordinateRotation::operator=(const dictionary& rhs)
202     if (debug)
203     {
204         Pout<< "coordinateRotation::operator=(const dictionary&) : "
205             << "assign from " << rhs << endl;
206     }
208     // allow as embedded sub-dictionary "coordinateRotation"
209     const dictionary& dict =
210     (
211         rhs.found(typeName_())
212       ? rhs.subDict(typeName_())
213       : rhs
214     );
216     vector axis1, axis2;
217     axisOrder order = e3e1;
219     if (dict.readIfPresent("e1", axis1) && dict.readIfPresent("e2", axis2))
220     {
221         order = e1e2;
222     }
223     else if
224     (
225         dict.readIfPresent("e2", axis1)
226      && dict.readIfPresent("e3", axis2)
227     )
228     {
229         order = e2e3;
230     }
231     else if
232     (
233         dict.readIfPresent("e3", axis1)
234      && dict.readIfPresent("e1", axis2)
235     )
236     {
237         order = e3e1;
238     }
239     else if (dict.found("axis") || dict.found("direction"))
240     {
241         // let it bomb if only one of axis/direction is defined
242         order = e3e1;
243         axis1 = vector(dict.lookup("axis"));
244         axis2 = vector(dict.lookup("direction"));
245     }
246     else
247     {
248         // unspecified axes revert to the global system
249         tensor::operator=(sphericalTensor::I);
250         return;
251     }
253     calcTransform(axis1, axis2, order);
257 // ************************************************************************* //