transformation3d.h

Go to the documentation of this file.
00001 #ifndef __Transformation3d_h
00002 #define __Transformation3d_h
00003 
00004 #include  "nmatrix.h"
00005 #include  "vertex.h"
00006 
00007 /*-----------------------------------------------------*/
00008 
00009 namespace Gascoigne
00010 {
00011 template<class BASE>
00012 class Transformation3d
00013 {
00014  protected:
00015 
00016   typedef nmatrix<double>    Matrix;
00017 
00018   BASE            B;
00019   mutable Matrix  X;
00020   mutable Matrix  dt, dti;
00021 
00022   inline void ComputeDT() const;
00023             // second derivatives tensor
00024   const nvector<Matrix>  ComputeDDT (const Vertex3d& xi) const;
00025 
00026  public:
00027 
00028   Transformation3d();
00029 
00030   const Matrix& DT () const {return dt ;}
00031   const Matrix& DTI() const {return dti;}
00032 
00033             // inverse of second derivatives tensor
00034   const nvector<Matrix>  DDTI(const Vertex3d& xi) const;
00035   
00036   inline double        J     () const;
00037   inline double        G     () const;
00038   inline Vertex3d      x     () const;
00039   inline Vertex3d      normal() const;
00040   inline void  init          (const Matrix& M) {X=M;}
00041   inline void  ReInit          (const Matrix& M) const {X=M;}
00042   inline void  point         (const Vertex3d& xi) const;
00043   inline void  point_boundary(int ie, const Vertex2d& s) const;
00044 };
00045 
00046 /*-----------------------------------------------------*/
00047 
00048 template<class BASE>
00049 inline Transformation3d<BASE>::Transformation3d() : B()
00050 {
00051   X  .memory(3,B.n());
00052   dt .memory(3,3);
00053   dti.memory(3,3);
00054 }
00055 
00056 /*-----------------------------------------------------*/
00057 
00058 template<class BASE>
00059 inline Vertex3d  Transformation3d<BASE>::x() const 
00060 {
00061   Vertex3d xp;
00062   for(int i=0;i<B.n();i++)
00063     {
00064       xp.x() += X(0,i) * B.phi(i);
00065       xp.y() += X(1,i) * B.phi(i);
00066       xp.z() += X(2,i) * B.phi(i);
00067     }
00068   return xp;
00069 }
00070 
00071 /*-----------------------------------------------------*/
00072 
00073 template<class BASE>
00074 inline Vertex3d  Transformation3d<BASE>::normal() const 
00075 {
00076   Vertex3d xn;
00077   dti.mult(xn,*B.normal3d());
00078   double xx = sqrt(xn*xn);
00079   xn /= xx;
00080   return xn;
00081 }
00082 
00083 /*-----------------------------------------------------*/
00084 
00085 template<class BASE>
00086 inline void  Transformation3d<BASE>::ComputeDT() const
00087 {
00088   dt.zero();
00089   for(int i=0;i<B.n();i++)
00090     {
00091       dt(0,0) += X(0,i) * B.phi_x(i);
00092       dt(0,1) += X(0,i) * B.phi_y(i);
00093       dt(0,2) += X(0,i) * B.phi_z(i);
00094       dt(1,0) += X(1,i) * B.phi_x(i);
00095       dt(1,1) += X(1,i) * B.phi_y(i);
00096       dt(1,2) += X(1,i) * B.phi_z(i);
00097       dt(2,0) += X(2,i) * B.phi_x(i);
00098       dt(2,1) += X(2,i) * B.phi_y(i);
00099       dt(2,2) += X(2,i) * B.phi_z(i);
00100     }
00101   dti(0,0) = dt(0,0);
00102   dti(0,1) = dt(1,0);
00103   dti(0,2) = dt(2,0);
00104   dti(1,0) = dt(0,1);
00105   dti(1,1) = dt(1,1);
00106   dti(1,2) = dt(2,1);
00107   dti(2,0) = dt(0,2);
00108   dti(2,1) = dt(1,2);
00109   dti(2,2) = dt(2,2);
00110 
00111   dti.gauss_jordan();
00112 }
00113 
00114 
00115 /*-----------------------------------------------------*/
00116 
00117 template<class BASE>
00118 inline const nvector<nmatrix<double> > Transformation3d<BASE>::DDTI(const Vertex3d& xi) const
00119 {
00120   const nvector<nmatrix<double> >& ddt = ComputeDDT(xi);
00121   
00122   nvector<nmatrix<double> > ddti(3,nmatrix<double> (3,3));
00123   Matrix dti_ = dti;
00124   
00125 
00126   dti_.transpose();
00127   for (int i=0;i<3;++i)
00128     {
00129       nmatrix<double> tmp(3,3);
00130 
00131       ddti[i].zero();
00132       for (int d=0;d<3;++d)
00133         ddti[i].add(-dti_(d,i),ddt[d]);
00134 
00135       ddti[i].mmult(tmp,dti_);
00136       dti_.mmult(ddti[i],tmp);
00137     }
00138   return ddti;  
00139 }
00140 
00141 /*-----------------------------------------------------*/
00142 template<class BASE>
00143 inline const nvector<nmatrix<double> > Transformation3d<BASE>::ComputeDDT(const Vertex3d& xi) const
00144 {
00145   const_cast<BASE*> (&B)->point(xi);
00146   
00147   nvector<nmatrix<double> > ddt(3,nmatrix<double> (3,3));
00148   for (int i=0;i<3;++i) ddt[i].zero();
00149 
00150   for (int i=0;i<B.n();++i)
00151     {
00152       for (int j=0;j<3;++j)
00153         {
00154           ddt[0](j,0) += X(j,i) * B.phi_xx(i);
00155           ddt[0](j,1) += X(j,i) * B.phi_xy(i);
00156           ddt[0](j,2) += X(j,i) * B.phi_xz(i);
00157 
00158           ddt[1](j,0) += X(j,i) * B.phi_xy(i);
00159           ddt[1](j,1) += X(j,i) * B.phi_yy(i);
00160           ddt[1](j,2) += X(j,i) * B.phi_yz(i);
00161           
00162           ddt[2](j,0) += X(j,i) * B.phi_xz(i);
00163           ddt[2](j,1) += X(j,i) * B.phi_yz(i);
00164           ddt[2](j,2) += X(j,i) * B.phi_zz(i);
00165         }
00166     }
00167   return ddt;  
00168 }
00169 
00170 /*-----------------------------------------------------*/
00171 
00172 template<class BASE>
00173 inline void  Transformation3d<BASE>::point(const Vertex3d& xi) const
00174 {
00175   B.point(xi);
00176   ComputeDT();
00177 }
00178 
00179 /*-----------------------------------------------------*/
00180 
00181 template<class BASE>
00182 inline void  Transformation3d<BASE>::point_boundary(int ie, const Vertex2d& s) const
00183 {
00184   B.point_boundary(ie,s);
00185   ComputeDT();
00186 }
00187 
00188 /*-----------------------------------------------------*/
00189 
00190 template<class BASE>
00191 inline double  Transformation3d<BASE>::J() const  
00192 {
00193   double h = dt.det();
00194   if (h<0)
00195     {
00196       std::cout << "Transformation3d<BASE>::J()  h = " << h << std::endl;
00197     }
00198   return h;
00199 }
00200 
00201 /*-----------------------------------------------------*/
00202 
00203 template<class BASE>
00204 inline double  Transformation3d<BASE>::G() const  
00205 {
00206   double d1phi=0,d2phi=0,d12phi=0;
00207   const fixarray<2,int>& fc = *B.faces();
00208   for (int i=0;i<3;++i)
00209     {
00210       d1phi+=dt(i,fc[0])*dt(i,fc[0]);
00211       d2phi+=dt(i,fc[1])*dt(i,fc[1]);
00212       d12phi+=dt(i,fc[0])*dt(i,fc[1]);
00213     }
00214   double h = d1phi*d2phi-d12phi*d12phi;
00215   assert(h>=0);
00216   return sqrt(h);
00217 }
00218 }
00219 
00220 #endif

Generated on Thu Sep 14 10:34:36 2006 for Gascoigne by  doxygen 1.4.7