transformation2d.h

Go to the documentation of this file.
00001 #ifndef __Transformation2d_h
00002 #define __Transformation2d_h
00003 
00004 #include  "nmatrix.h"
00005 #include  "vertex.h"
00006 
00007 /*-----------------------------------------------------*/
00008 
00009 namespace Gascoigne
00010 {
00011 template<class BASE>
00012 class Transformation2d
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 
00024             // second derivatives tensor
00025   const nvector<Matrix>  ComputeDDT (const Vertex2d& xi) const;
00026 
00027  public:
00028 
00029   Transformation2d();
00030 
00031   const Matrix& DT () const {return dt ;}
00032   const Matrix& DTI() const {return dti;}
00033 
00034             // inverse of second derivatives tensor
00035   const nvector<Matrix>  DDTI(const Vertex2d& xi) const;
00036   
00037   inline double        J     () const;
00038   inline double        G     () const;
00039   inline Vertex2d      x     () const;
00040   inline Vertex2d      normal() const;
00041   inline void  init          (const Matrix& M) {X=M;}
00042   inline void  ReInit          (const Matrix& M) const {X=M;}
00043   inline void  point         (const Vertex2d& xi) const;
00044   inline void  point_boundary(int ie, const Vertex1d& s) const;
00045 };
00046 
00047 /*-----------------------------------------------------*/
00048 
00049 template<class BASE>
00050 inline Transformation2d<BASE>::Transformation2d() : B()
00051 {
00052   int n = B.n();
00053   X.memory(2,n);
00054   dt .memory(2,2);
00055   dti.memory(2,2);
00056 }
00057 
00058 /*-----------------------------------------------------*/
00059 
00060 template<class BASE>
00061 inline Vertex2d  Transformation2d<BASE>::x() const 
00062 {
00063   Vertex2d xp;
00064   for(int i=0;i<B.n();i++)
00065     {
00066       xp.x() += X(0,i) * B.phi(i);
00067       xp.y() += X(1,i) * B.phi(i);
00068     }
00069   return xp;
00070 }
00071 
00072 /*-----------------------------------------------------*/
00073 
00074 template<class BASE>
00075 inline Vertex2d  Transformation2d<BASE>::normal() const 
00076 {
00077   Vertex2d xn;
00078   dti.mult(xn,*B.normal2d());
00079   double xx = sqrt(xn*xn);
00080   xn /= xx;
00081   return xn;
00082 }
00083 
00084 /*-----------------------------------------------------*/
00085 
00086 template<class BASE>
00087 inline void  Transformation2d<BASE>::ComputeDT() const
00088 {
00089   dt.zero();
00090   for(int i=0;i<B.n();i++)
00091     {
00092       dt(0,0) += X(0,i) * B.phi_x(i);
00093       dt(0,1) += X(0,i) * B.phi_y(i);
00094       dt(1,0) += X(1,i) * B.phi_x(i);
00095       dt(1,1) += X(1,i) * B.phi_y(i);
00096     }
00097   dti(0,0) = dt(0,0);      
00098   dti(0,1) = dt(1,0);      
00099   dti(1,0) = dt(0,1);      
00100   dti(1,1) = dt(1,1);      
00101   dti.gauss_jordan();
00102 }
00103 
00104 /*-----------------------------------------------------*/
00105 
00106 template<class BASE>
00107 inline const nvector<nmatrix<double> > Transformation2d<BASE>::ComputeDDT(const Vertex2d& xi) const
00108 {
00109   const_cast<BASE*> (&B)->point(xi);
00110   
00111   nvector<nmatrix<double> > ddt(2,nmatrix<double> (2,2));
00112   for (int i=0;i<2;++i) ddt[i].zero();
00113 
00114   for (int i=0;i<B.n();++i)
00115     {
00116       ddt[0](0,0) += X(0,i) * B.phi_xx(i);
00117       ddt[0](0,1) += X(0,i) * B.phi_xy(i);
00118       ddt[0](1,0) += X(1,i) * B.phi_xx(i);
00119       ddt[0](1,1) += X(1,i) * B.phi_xy(i);
00120 
00121       ddt[1](0,0) += X(0,i) * B.phi_xy(i);
00122       ddt[1](0,1) += X(0,i) * B.phi_yy(i);
00123       ddt[1](1,0) += X(1,i) * B.phi_xy(i);
00124       ddt[1](1,1) += X(1,i) * B.phi_yy(i);
00125     }
00126   return ddt;  
00127 }
00128 
00129 /*-----------------------------------------------------*/
00130 
00131 template<class BASE>
00132 inline const nvector<nmatrix<double> > Transformation2d<BASE>::DDTI(const Vertex2d& xi) const
00133 {
00134   const nvector<nmatrix<double> >& ddt = ComputeDDT(xi);
00135   
00136   nvector<nmatrix<double> > ddti(2,nmatrix<double> (2,2));
00137   Matrix dti_ = dti;
00138   
00139 
00140   dti_.transpose();
00141   for (int i=0;i<2;++i)
00142     {
00143       nmatrix<double> tmp(2,2);
00144 
00145       ddti[i].zero();
00146       for (int d=0;d<2;++d)
00147         ddti[i].add(-dti_(d,i),ddt[d]);
00148       
00149       ddti[i].mmult(tmp,dti_);
00150       dti_.mmult(ddti[i],tmp);
00151     }
00152   return ddti;  
00153 }
00154 
00155 /*-----------------------------------------------------*/
00156 
00157 template<class BASE>
00158 inline void  Transformation2d<BASE>::point(const Vertex2d& xi) const
00159 {
00160   B.point(xi);
00161   ComputeDT();
00162 }
00163 
00164 /*-----------------------------------------------------*/
00165 
00166 template<class BASE>
00167 inline void  Transformation2d<BASE>::point_boundary(int ie, const Vertex1d& s) const
00168 {
00169   B.point_boundary(ie,s);
00170   ComputeDT();
00171 }
00172 
00173 /*-----------------------------------------------------*/
00174 
00175 template<class BASE>
00176 inline double  Transformation2d<BASE>::J() const  
00177 {
00178   return dt(0,0)*dt(1,1)-dt(1,0)*dt(0,1);
00179 }
00180 
00181 /*-----------------------------------------------------*/
00182 
00183 template<class BASE>
00184 inline double  Transformation2d<BASE>::G() const  
00185 {
00186   Vertex2d xt;
00187   dt.mult(xt,*B.tangent2d());
00188   return xt.norm_l2();
00189 }
00190 }
00191 
00192 #endif

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