FreeFEM3D (aka ff3d): ConformTransformation.hpp Source File

00001 //  This file is part of ff3d - http://www.freefem.org/ff3d
00002 //  Copyright (C) 2001, 2002, 2003 St�phane Del Pino
00003 
00004 //  This program is free software; you can redistribute it and/or modify
00005 //  it under the terms of the GNU General Public License as published by
00006 //  the Free Software Foundation; either version 2, or (at your option)
00007 //  any later version.
00008 
00009 //  This program is distributed in the hope that it will be useful,
00010 //  but WITHOUT ANY WARRANTY; without even the implied warranty of
00011 //  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00012 //  GNU General Public License for more details.
00013 
00014 //  You should have received a copy of the GNU General Public License
00015 //  along with this program; if not, write to the Free Software Foundation,
00016 //  Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.  
00017 
00018 //  $Id: ConformTransformation.hpp,v 1.18 2007/08/29 08:37:14 yakoubix Exp $
00019 
00020 #ifndef CONFORM_TRANSFORMATION_HPP
00021 #define CONFORM_TRANSFORMATION_HPP
00022 
00023 #include <TinyVector.hpp>
00024 #include <TinyMatrix.hpp>
00025 
00026 #include <Vertex.hpp>
00027 #include <Cell.hpp>
00028 
00029 #include <Hexahedron.hpp>
00030 #include <CartesianHexahedron.hpp>
00031 
00032 #include <Tetrahedron.hpp>
00033 
00034 #include <Triangle.hpp>
00035 #include <Quadrangle.hpp>
00036 
00037 #include <QuadratureFormula.hpp>
00038 
00039 #include <ErrorHandler.hpp>
00040 
00041 class Domain;
00042 class ScalarFunctionBase;
00043 
00044 class ConformTransformationP1Triangle
00045 {
00046 public:
00047   friend class ConformTransformationP1TriangleJacobian;
00048 
00049 private:
00050   const Triangle& __T;
00051 
00052   const TinyVector<3, real_t> __a;
00053   const TinyVector<3, real_t> __b_a; // (b-a) vector
00054   const TinyVector<3, real_t> __c_a; // (c-a) vector
00055 
00056 public:
00058   inline void value(const TinyVector<3, real_t>& X,
00059                     TinyVector<3, real_t>& result) const
00060   {
00061     value(X[0],X[1],X[2],result);
00062   }
00063 
00065   inline void value(const real_t& x, const real_t& y, const real_t& z,
00066                     TinyVector<3, real_t>& result) const
00067   {
00068     result  = __a;
00069     result += __b_a*x;
00070     result += __c_a*y;
00071   }
00072 
00073 
00074 
00075   inline void dx(const TinyVector<3, real_t>& X,
00076                  TinyVector<3, real_t>& __result) const
00077   {
00078     dx(X[0], X[1], X[2], __result);
00079   }
00080 
00081   inline void dx(const real_t& x, const real_t& y, const real_t& z,
00082                  TinyVector<3, real_t>& __result) const
00083   {
00084     __result = __b_a;
00085   }
00086 
00087   inline void dy(const TinyVector<3, real_t>& X,
00088                  TinyVector<3, real_t>& __result) const
00089   {
00090     dy(X[0], X[1], X[2], __result);
00091   }
00092 
00093   inline void dy(const real_t& x, const real_t& y, const real_t& z,
00094                  TinyVector<3, real_t>& __result) const
00095   {
00096     __result = __c_a;
00097   }
00098 
00099   inline void dz(const TinyVector<3, real_t>& X,
00100                  TinyVector<3, real_t>& __result) const
00101   {
00102     dz(X[0], X[1], X[2], __result);
00103   }
00104 
00105   inline void dz(const real_t& x, const real_t& y, const real_t& z,
00106                  TinyVector<3, real_t>& __result) const
00107   {
00108     __result = 0;
00109   }
00110 
00111 
00112 
00113 
00114 
00115 
00117   inline bool invertT(const TinyVector<3, real_t>& X,
00118                       TinyVector<3, real_t>& Xhat) const
00119   {
00120     return invertT(X[0],X[1],X[2], Xhat);
00121   }
00122 
00124 
00136   inline bool invertT(const real_t& x, const real_t& y, const real_t& z,
00137                       TinyVector<3, real_t>& Xhat) const
00138   {
00139     TinyMatrix<2,2,real_t> A;
00140     for (size_t i=0; i<2; ++i) {
00141       A(i,0) = __b_a[i];
00142       A(i,1) = __c_a[i];
00143     }
00144 
00145     TinyVector<2,real_t> b;
00146     b[0] = x-__a[0];
00147     b[1] = y-__a[1];
00148 
00149     TinyVector<2,real_t> u = b/A;
00150     Xhat[0] = u[0];
00151     Xhat[1] = u[1];
00152     Xhat[2] = 0;
00153 
00154     TinyVector<3,real_t> result  = __a;
00155     result += __b_a*u[0];
00156     result += __c_a*u[1];
00157 
00158     return true;
00159   }
00160 
00164   real_t integrate(const ScalarFunctionBase& f) const;
00165 
00166   ConformTransformationP1Triangle(const Triangle& T)
00167     : __T(T),
00168       __a(T(0)),
00169       __b_a(T(1)-T(0)),
00170       __c_a(T(2)-T(0))
00171   {
00172     ;
00173   }
00174 };
00175 
00176 
00177 class ConformTransformationP1TriangleJacobian
00178 {
00179 public:
00180   typedef ConformTransformationP1Triangle AssociatedTransformation;
00181 
00182 private:
00183   const ConformTransformationP1Triangle& __T;
00184 
00185   TinyMatrix<3,3> __invJacobian;
00186 
00187   real_t __det;
00188 
00189   void dx(const real_t& x, const real_t& y, const real_t& z,
00190           TinyVector<3, real_t>& __result) const
00191   {
00192     //     __result = __T.__h[0]*x;
00193   }
00194 
00195   void dy(const real_t& x, const real_t& y, const real_t& z,
00196           TinyVector<3, real_t>& __result) const
00197   {
00198     //     __result = __T.__h[1]*y;
00199   }
00200 
00201   void dz(const real_t& x, const real_t& y, const real_t& z,
00202           TinyVector<3, real_t>& __result) const
00203   {
00204     //     __result = __T.__h[2]*z;
00205   }
00206   
00207   
00208 public:
00209   
00213   const TinyMatrix<3,3,real_t>& invJacobian() const
00214   {
00215     return __invJacobian;
00216   }
00217 
00219   const real_t& invJacobian(const size_t i, const size_t j) const
00220   {
00221     return __invJacobian(i,j);
00222   }
00223 
00225   const real_t& jacobianDet() const
00226   {
00227     return __det;
00228   }
00229 
00230   ConformTransformationP1TriangleJacobian(const ConformTransformationP1Triangle& T)
00231     : __T(T),
00232       __invJacobian(0),
00233       __det(0)
00234   {
00235 //     throw ErrorHandler(__FILE__,__LINE__,
00236 //                     "not implemented",
00237 //                     ErrorHandler::unexpected);
00238 
00239     //     for (size_t i=0; i<3; ++i)
00240     //       __invJacobian(i,i) = 1./(__T.__h[i]);
00241 
00242     __det = 2*__T.__T.volume();
00243   }
00244 
00245 };
00246 
00247 class ConformTransformationQ1Quadrangle
00248 {
00249 public:
00250   friend class ConformTransformationQ1QuadrangleJacobian;
00251 
00252 private:
00253   const Quadrangle& __Q;
00254 
00255   const TinyVector<3, real_t> __a;
00256   const TinyVector<3, real_t> __b_a;    // (b-a) vector
00257   const TinyVector<3, real_t> __d_a;    // (d-a) vector
00258   const TinyVector<3, real_t> __ca_b_d; // (c+a-b-d) vector
00259 
00260 public:
00262   inline void value(const TinyVector<3, real_t>& X,
00263                     TinyVector<3, real_t>& result) const
00264   {
00265     value(X[0],X[1],X[2],result);
00266   }
00267 
00269   inline void value(const real_t& x, const real_t& y, const real_t& z,
00270                     TinyVector<3, real_t>& result) const
00271   {
00272     result  = __a;
00273     result += __b_a*x;
00274     result += __d_a*y;
00275     result += __ca_b_d*x*y;
00276   }
00277 
00279   inline bool invertT(const TinyVector<3, real_t>& X,
00280                       TinyVector<3, real_t>& Xhat) const
00281   {
00282     return invertT(X[0],X[1],X[2], Xhat);
00283   }
00284 
00286   inline bool invertT(const real_t& x, const real_t& y, const real_t& z,
00287                       TinyVector<3, real_t>& Xhat) const
00288   {
00289     // Use least square to solve the linear system
00290     throw ErrorHandler(__FILE__,__LINE__,
00291                        "not implemented",
00292                        ErrorHandler::unexpected);
00293 
00294     return false;
00295   }
00296 
00300   real_t integrate(const ScalarFunctionBase& f) const;
00301 
00302   ConformTransformationQ1Quadrangle(const Quadrangle& Q)
00303     : __Q(Q),
00304       __a(Q(0)),
00305       __b_a(Q(1)-Q(0)),
00306       __d_a(Q(3)-Q(0)),
00307       __ca_b_d(Q(2)-__b_a-Q(3))
00308   {
00309     ;
00310   }
00311 };
00312 
00313 class ConformTransformationQ1QuadrangleJacobian
00314 {
00315 public:
00316   typedef ConformTransformationQ1Quadrangle AssociatedTransformation;
00317 
00318 private:
00319   const ConformTransformationQ1Quadrangle& __T;
00320 
00321   TinyMatrix<3,3> __invJacobian;
00322 
00323   real_t __det;
00324 
00325   void dx(const real_t& x, const real_t& y, const real_t& z,
00326           TinyVector<3, real_t>& __result) const
00327   {
00328     //     __result = __T.__h[0]*x;
00329   }
00330 
00331   void dy(const real_t& x, const real_t& y, const real_t& z,
00332           TinyVector<3, real_t>& __result) const
00333   {
00334     //     __result = __T.__h[1]*y;
00335   }
00336 
00337   void dz(const real_t& x, const real_t& y, const real_t& z,
00338           TinyVector<3, real_t>& __result) const
00339   {
00340     //     __result = __T.__h[2]*z;
00341   }
00342 
00343 
00344 public:
00345 
00349   const TinyMatrix<3,3,real_t>& invJacobian() const
00350   {
00351     return __invJacobian;
00352   }
00353 
00355   const real_t& invJacobian(const size_t i, const size_t j) const
00356   {
00357     return __invJacobian(i,j);
00358   }
00359 
00361   const real_t& jacobianDet() const
00362   {
00363     return __det;
00364   }
00365 
00366   ConformTransformationQ1QuadrangleJacobian(const ConformTransformationQ1Quadrangle& T)
00367     : __T(T),
00368       __invJacobian(0),
00369       __det(0)
00370   {
00371     __det = __T.__Q.volume();
00372   }
00373 };
00374 
00375 class ConformTransformationQ1Hexahedron
00376 {
00377 public:
00378   friend class ConformTransformationQ1HexahedronJacobian;
00379   typedef ConformTransformationQ1Quadrangle
00380   BoundaryConformTransformation;
00381 
00382 private:
00383   const Hexahedron& __H;
00384 
00385   TinyVector<3, real_t>& __a;
00386   TinyVector<3, real_t>& __b;
00387   TinyVector<3, real_t>& __c;
00388   TinyVector<3, real_t>& __d;
00389   TinyVector<3, real_t>& __e;
00390   TinyVector<3, real_t>& __f;
00391   TinyVector<3, real_t>& __g;
00392   TinyVector<3, real_t>& __h;
00393 
00394   TinyVector<8, TinyVector<3, real_t> > __vertices;
00395 
00396   bool __inside(const real_t& x, const real_t& y, const real_t& z,
00397                 TinyVector<6, bool>& faces) const;
00398 
00399 public:
00400 
00401   inline void dx(const TinyVector<3, real_t>& X,
00402                  TinyVector<3, real_t>& __result) const
00403   {
00404     dx(X[0], X[1], X[2], __result);
00405   }
00406 
00407   inline void dx(const real_t& x, const real_t& y, const real_t& z,
00408                  TinyVector<3, real_t>& __result) const
00409   {
00410     __result = __b;
00411     __result += y*__e;
00412     __result += z*__f;
00413     __result += (y*z)*__h;
00414   }
00415 
00416   inline void dy(const TinyVector<3, real_t>& X,
00417                  TinyVector<3, real_t>& __result) const
00418   {
00419     dy(X[0], X[1], X[2], __result);
00420   }
00421 
00422   inline void dy(const real_t& x, const real_t& y, const real_t& z,
00423                  TinyVector<3, real_t>& __result) const
00424   {
00425     __result = __c;
00426     __result += x*__e;
00427     __result += z*__g;
00428     __result += (x*z)*__h;
00429   }
00430 
00431   inline void dz(const TinyVector<3, real_t>& X,
00432                  TinyVector<3, real_t>& __result) const
00433   {
00434     dz(X[0], X[1], X[2], __result);
00435   }
00436 
00437   inline void dz(const real_t& x, const real_t& y, const real_t& z,
00438                  TinyVector<3, real_t>& __result) const
00439   {
00440     __result = __d;
00441     __result += x*__f;
00442     __result += y*__g;
00443     __result += (x*y)*__h;
00444   }
00445 
00446   inline void value(const TinyVector<3, real_t>& X,
00447                     TinyVector<3, real_t>& __result) const
00448   {
00449     value(X[0],X[1],X[2],__result);
00450   }
00451 
00452   inline void value(const real_t& x, const real_t& y, const real_t& z,
00453                     TinyVector<3, real_t>& __result) const
00454   {
00455     __result  = __a;
00456     __result += x*__b;
00457     __result += y*__c;
00458     __result += z*__d;
00459     __result += (x*y)*__e;
00460     __result += (x*z)*__f;
00461     __result += (y*z)*__g;
00462     __result += (x*y*z)*__h;
00463   }
00464 
00468   bool invertT(const real_t& x, const real_t& y, const real_t& z,
00469                TinyVector<3, real_t>& Xhat) const;
00470 
00472   inline bool invertT(const TinyVector<3, real_t>& X,
00473                       TinyVector<3, real_t>& Xhat) const
00474   {
00475     return invertT(X[0],X[1],X[2], Xhat);
00476   }
00477 
00482   real_t integrate(const ScalarFunctionBase& f) const;
00483 
00484   ConformTransformationQ1Hexahedron(const Hexahedron& H)
00485     : __H(H),
00486       __a(__vertices[0]),
00487       __b(__vertices[1]),
00488       __c(__vertices[2]),
00489       __d(__vertices[3]),
00490       __e(__vertices[4]),
00491       __f(__vertices[5]),
00492       __g(__vertices[6]),
00493       __h(__vertices[7])
00494   {
00495     TinyVector<8,TinyVector<3, real_t> > jCoefs;
00496 
00497     jCoefs[0]=H(0);
00498     jCoefs[1]=H(1) - jCoefs[0];
00499     jCoefs[2]=H(3) - jCoefs[0];
00500     jCoefs[3]=H(4) - jCoefs[0];
00501     jCoefs[4]=H(2) + jCoefs[0] - H(1) - H(3);
00502     jCoefs[5]=H(5) + jCoefs[0] - H(1) - H(4);
00503     jCoefs[6]=H(7) + jCoefs[0] - H(3) - H(4);
00504     jCoefs[7]
00505       = H(1) - jCoefs[0] - H(2) + H(3)
00506       + H(4) - H(5)
00507       + H(6) - H(7);
00508 
00509     __vertices = jCoefs;
00510   }
00511 };
00512 
00513 class ConformTransformationQ1HexahedronJacobian
00514 {
00515 private:
00516   const ConformTransformationQ1Hexahedron& __T;
00517 
00518   TinyMatrix<3,3> __jacobian;
00519   TinyMatrix<3,3> __invJacobian;
00520 
00521   real_t __det;
00522 
00523 public:
00524 
00526   const TinyMatrix<3,3,real_t>& jacobian() const
00527   {
00528     return __jacobian;
00529   }
00530 
00534   const TinyMatrix<3,3,real_t>& invJacobian() const
00535   {
00536     return __invJacobian;
00537   }
00538 
00540   const real_t& jacobianDet() const
00541   {
00542     return __det;
00543   }
00544 
00546   const real_t& invJacobian(const size_t i, const size_t j) const
00547   {
00548     return __invJacobian(i,j);
00549   }
00550 
00551   ConformTransformationQ1HexahedronJacobian(const ConformTransformationQ1Hexahedron& T)
00552     : __T(T)
00553   {
00554     const TinyVector<QuadratureFormulaQ1Hexahedron::numberOfQuadraturePoints,
00555                      TinyVector<3, real_t> >& IntegrationVertices
00556       = QuadratureFormulaQ1Hexahedron::instance().vertices();
00557 
00558     __jacobian = 0;
00559 
00560     TinyVector<QuadratureFormulaQ1Hexahedron::numberOfQuadraturePoints, TinyVector<3, real_t> > dxPhi;
00561     TinyVector<QuadratureFormulaQ1Hexahedron::numberOfQuadraturePoints, TinyVector<3, real_t> > dyPhi;
00562     TinyVector<QuadratureFormulaQ1Hexahedron::numberOfQuadraturePoints, TinyVector<3, real_t> > dzPhi;
00563 
00564     for (size_t i=0; i<QuadratureFormulaQ1Hexahedron::numberOfQuadraturePoints; ++i) {
00565       __T.dx(IntegrationVertices[i][0],
00566              IntegrationVertices[i][1],
00567              IntegrationVertices[i][2],dxPhi[i]);
00568       __T.dy(IntegrationVertices[i][0],
00569              IntegrationVertices[i][1],
00570              IntegrationVertices[i][2],dyPhi[i]);
00571       __T.dz(IntegrationVertices[i][0],
00572              IntegrationVertices[i][1],
00573              IntegrationVertices[i][2],dzPhi[i]);
00574       for (size_t j=0; j<3; ++j) {
00575         __jacobian(0,j) += dxPhi[i][j];
00576         __jacobian(1,j) += dyPhi[i][j];
00577         __jacobian(2,j) += dzPhi[i][j];
00578       }
00579     }
00580 
00581     __jacobian *= 1./QuadratureFormulaQ1Hexahedron::numberOfQuadraturePoints;
00582 
00583     TinyVector<3, TinyVector<3, real_t> > I;
00584     for (size_t i=0; i<3; ++i) {
00585       for (size_t j=0; j<3; ++j)
00586         I[i][j] = (i==j)?1:0;
00587     }
00588 
00589     TinyVector<3,  TinyVector<3, real_t> > J;
00590     __det = ::gaussPivot(__jacobian,I,J);
00591 
00592     for (size_t i=0; i<3; ++i) {
00593       for (size_t j=0; j<3; ++j) {
00595         __invJacobian(j,i) = J[i][j];
00596       }
00597     }
00598   }
00599 };
00600 
00601 
00602 class ConformTransformationQ1CartesianHexahedron
00603 {
00604 public:
00605   friend class ConformTransformationQ1CartesianHexahedronJacobian;
00606 
00607   typedef ConformTransformationQ1Quadrangle BoundaryConformTransformation;
00608 
00609 private:
00610   const CartesianHexahedron& __CH;
00611 
00612   const TinyVector<3, real_t> __a;
00613   const TinyVector<3, real_t> __h;
00614 
00615 public:
00616   inline void dx(const TinyVector<3, real_t>& X,
00617                  TinyVector<3, real_t>& __result) const
00618   {
00619     dx(X[0], X[1], X[2], __result);
00620   }
00621 
00622   inline void dx(const real_t& x, const real_t& y, const real_t& z,
00623                  TinyVector<3, real_t>& __result) const
00624   {
00625     __result[0] = __h[0];
00626     __result[1] = 0;
00627     __result[2] = 0;
00628   }
00629 
00630   inline void dy(const TinyVector<3, real_t>& X,
00631                  TinyVector<3, real_t>& __result) const
00632   {
00633     dy(X[0], X[1], X[2], __result);
00634   }
00635 
00636   inline void dy(const real_t& x, const real_t& y, const real_t& z,
00637                  TinyVector<3, real_t>& __result) const
00638   {
00639     __result[0] = 0;
00640     __result[1] = __h[1];
00641     __result[2] = 0;
00642   }
00643 
00644   inline void dz(const TinyVector<3, real_t>& X,
00645                  TinyVector<3, real_t>& __result) const
00646   {
00647     dz(X[0], X[1], X[2], __result);
00648   }
00649 
00650   inline void dz(const real_t& x, const real_t& y, const real_t& z,
00651                  TinyVector<3, real_t>& __result) const
00652   {
00653     __result[0] = 0;
00654     __result[1] = 0;
00655     __result[2] = __h[2];
00656   }
00657 
00659   inline void value(const TinyVector<3, real_t>& X,
00660                     TinyVector<3, real_t>& result) const
00661   {
00662     value(X[0],X[1],X[2],result);
00663   }
00664 
00666   inline void value(const real_t& x, const real_t& y, const real_t& z,
00667                     TinyVector<3, real_t>& result) const
00668   {
00669     result  = __a;
00670     result[0] += __h[0]*x;
00671     result[1] += __h[1]*y;
00672     result[2] += __h[2]*z;
00673   }
00674 
00676   inline bool invertT(const TinyVector<3, real_t>& X,
00677                       TinyVector<3, real_t>& Xhat) const
00678   {
00679     return invertT(X[0],X[1],X[2], Xhat);
00680   }
00681 
00683   inline bool invertT(const real_t& x, const real_t& y, const real_t& z,
00684                       TinyVector<3, real_t>& Xhat) const
00685   {
00686     Xhat[0] = (x-__a[0])/__h[0];
00687     Xhat[1] = (y-__a[1])/__h[1];
00688     Xhat[2] = (z-__a[2])/__h[2];
00689     return true;
00690   }
00691 
00695   real_t integrate(const ScalarFunctionBase& f) const;
00696 
00701   real_t integrateCharacteristic(const Domain& d) const;
00702 
00703   ConformTransformationQ1CartesianHexahedron(const CartesianHexahedron& H)
00704     : __CH(H),
00705       __a(H(0)),
00706       __h(H(6)-H(0))
00707   {
00708     ;
00709   }
00710 };
00711 
00712 
00713 class ConformTransformationQ1CartesianHexahedronJacobian
00714 {
00715 public:
00716   typedef ConformTransformationQ1CartesianHexahedron AssociatedTransformation;
00717 
00718 private:
00719   ConformTransformationQ1CartesianHexahedron& __T;
00720 
00721   TinyMatrix<3,3> __jacobian;
00722 
00723   TinyMatrix<3,3> __invJacobian;
00724 
00725   real_t __det;
00726 
00727   void dx(const real_t& x, const real_t& y, const real_t& z,
00728           TinyVector<3, real_t>& __result) const
00729   {
00730     __result = __T.__h[0]*x;
00731   }
00732 
00733   void dy(const real_t& x, const real_t& y, const real_t& z,
00734           TinyVector<3, real_t>& __result) const
00735   {
00736     __result = __T.__h[1]*y;
00737   }
00738 
00739   void dz(const real_t& x, const real_t& y, const real_t& z,
00740           TinyVector<3, real_t>& __result) const
00741   {
00742     __result = __T.__h[2]*z;
00743   }
00744 public:
00745 
00746   inline real_t
00747   jacobian(const size_t i, const size_t j) const
00748   {
00749     return __jacobian(i,j);
00750   }
00751 
00755   const TinyMatrix<3,3,real_t>& invJacobian() const
00756   {
00757     return __invJacobian;
00758   }
00759 
00761   const real_t& invJacobian(const size_t i, const size_t j) const
00762   {
00763     return __invJacobian(i,j);
00764   }
00765 
00767   const real_t& jacobianDet() const
00768   {
00769     return __det;
00770   }
00771 
00772   ConformTransformationQ1CartesianHexahedronJacobian(ConformTransformationQ1CartesianHexahedron& T)
00773     : __T(T),
00774       __jacobian(0),
00775       __invJacobian(0),
00776       __det(0)
00777   {
00778     for (size_t i=0; i<3; ++i) {
00779       __invJacobian(i,i) = 1./(__T.__h[i]);
00780       __jacobian(i,i) = (__T.__h[i]);
00781     }
00782 
00783     __det = __T.__h[0] * __T.__h[1] * __T.__h[2];
00784 
00785   }
00786 };
00787 
00788 class ConformTransformationP1Tetrahedron
00789 {
00790 public:
00791   friend class ConformTransformationP1TetrahedronJacobian;
00792   typedef ConformTransformationP1Tetrahedron BoundaryConformTransformation;
00793 
00794 private:
00795   const Tetrahedron& __T;
00796 
00797   const TinyVector<3, real_t> __b;
00798   TinyMatrix<3,3, real_t> __A;
00799   TinyMatrix<3,3, real_t> __A_1;
00800 
00801 public:
00802   inline void dx(const TinyVector<3, real_t>& X,
00803                  TinyVector<3, real_t>& __result) const
00804   {
00805     dx(X[0], X[1], X[2], __result);
00806   }
00807 
00808   inline void dx(const real_t& x, const real_t& y, const real_t& z,
00809                  TinyVector<3, real_t>& __result) const
00810   {
00811     for (size_t i=0; i<3; ++i) __result[i] = __A(i,0);
00812   }
00813 
00814   inline void dy(const TinyVector<3, real_t>& X,
00815                  TinyVector<3, real_t>& __result) const
00816   {
00817     dy(X[0], X[1], X[2], __result);
00818   }
00819 
00820   inline void dy(const real_t& x, const real_t& y, const real_t& z,
00821                  TinyVector<3, real_t>& __result) const
00822   {
00823     for (size_t i=0; i<3; ++i) __result[i] = __A(i,1);
00824   }
00825 
00826   inline void dz(const TinyVector<3, real_t>& X,
00827                  TinyVector<3, real_t>& __result) const
00828   {
00829     dz(X[0], X[1], X[2], __result);
00830   }
00831 
00832   inline void dz(const real_t& x, const real_t& y, const real_t& z,
00833                  TinyVector<3, real_t>& __result) const
00834   {
00835     for (size_t i=0; i<3; ++i) __result[i] = __A(i,2);
00836   }
00837 
00839   inline void value(const TinyVector<3, real_t>& X,
00840                     TinyVector<3, real_t>& result) const
00841   {
00842     result = __A*X;
00843     result += __b;
00844   }
00845 
00847   inline void value(const real_t& x, const real_t& y, const real_t& z,
00848                     TinyVector<3, real_t>& result) const
00849   {
00850     TinyVector<3,real_t> X(x,y,z);
00851     value(X, result);
00852   }
00853 
00855   inline bool invertT(const TinyVector<3, real_t>& X,
00856                       TinyVector<3, real_t>& Xhat) const
00857   {
00858     Xhat = __A_1*(X-__b);
00859     return true;
00860   }
00861 
00863   inline bool invertT(const real_t& x, const real_t& y, const real_t& z,
00864                       TinyVector<3, real_t>& Xhat) const
00865   {
00866     TinyVector<3,real_t> X(x,y,z);
00867     return invertT(X, Xhat);
00868   }
00869 
00873   real_t integrate(const ScalarFunctionBase& f) const;
00874 
00875   ConformTransformationP1Tetrahedron(const Tetrahedron& T)
00876     : __T(T),
00877       __b(T(0))
00878   {
00879     for (size_t j=0; j<3; ++j) {
00880       for (size_t i=0; i<3; ++i) {
00881         __A(i,j) = T(j+1)[i]-T(0)[i];
00882       }
00883     }
00884 
00885     __A_1 = __A.invert();
00886   }
00887 };
00888 
00889 class ConformTransformationP1TetrahedronJacobian
00890 {
00891 public:
00892   typedef ConformTransformationP1Tetrahedron AssociatedTransformation;
00893 
00894 private:
00895   ConformTransformationP1Tetrahedron& __T;
00896 
00897   TinyMatrix<3,3>& __invJacobian;
00898 
00899   real_t __det;
00900 public:
00901 
00905   const TinyMatrix<3,3,real_t>& invJacobian() const
00906   {
00907     return __invJacobian;
00908   }
00909 
00911   const real_t& invJacobian(const size_t i, const size_t j) const
00912   {
00913     return __invJacobian(i,j);
00914   }
00915 
00917   const real_t& jacobianDet() const
00918   {
00919     return __det;
00920   }
00921 
00922   ConformTransformationP1TetrahedronJacobian(ConformTransformationP1Tetrahedron& T)
00923     : __T(T),
00924       __invJacobian(T.__A_1),
00925       __det(0)
00926   {
00927     __det = det(T.__A);
00928   }
00929 
00930 };
00931 
00932 #endif // CONFORM_TRANSFORMATION_HPP