dune-localfunctions  2.2.0
raviartthomas0q3dall.hh
Go to the documentation of this file.
00001 #ifndef DUNE_RT0Q3DALL_HH
00002 #define DUNE_RT0Q3DALL_HH
00003 
00004 #include <cstddef>
00005 #include <vector>
00006 
00007 #include <dune/common/fmatrix.hh>
00008 
00009 #include <dune/localfunctions/common/localbasis.hh>
00010 #include <dune/localfunctions/common/localkey.hh>
00011 
00012 namespace Dune 
00013 {
00022   template<class D, class R>
00023   class RT0Q3DLocalBasis
00024   {
00025   public:
00026         typedef LocalBasisTraits<D,3,Dune::FieldVector<D,3>,R,3,Dune::FieldVector<R,3>,
00027                                                            Dune::FieldMatrix<R,3,3> > Traits;
00028 
00030         RT0Q3DLocalBasis ()
00031         {
00032           sign0 = sign1 = sign2 = sign3 = sign4 = sign5 = 1.0;
00033         }
00034 
00036         RT0Q3DLocalBasis (unsigned int s)
00037         {
00038           sign0 = sign1 = sign2 = sign3 = sign4 = sign5 = 1.0;
00039           if (s&1) sign0 = -1.0;
00040           if (s&2) sign1 = -1.0;
00041           if (s&4) sign2 = -1.0;
00042           if (s&8) sign3 = -1.0;
00043           if (s&16) sign4 = -1.0;
00044           if (s&32) sign5 = -1.0;
00045         }
00046 
00048         unsigned int size () const
00049         {
00050           return 6;
00051         }
00052 
00054         inline void evaluateFunction (const typename Traits::DomainType& in,
00055                                                                   std::vector<typename Traits::RangeType>& out) const
00056         { 
00057           out.resize(6);
00058           out[0][0] = sign0*(in[0]-1.0); out[0][1]=0.0;               out[0][2]=0.0;
00059           out[1][0] = sign1*(in[0]);     out[1][1]=0.0;               out[1][2]=0.0;
00060           out[2][0] = 0.0;               out[2][1]=sign2*(in[1]-1.0); out[2][2]=0.0;
00061           out[3][0] = 0.0;               out[3][1]=sign3*(in[1]);     out[3][2]=0.0;
00062           out[4][0] = 0.0;               out[4][1]=0.0;               out[4][2]=sign4*(in[2]-1.0);
00063           out[5][0] = 0.0;               out[5][1]=0.0;               out[5][2]=sign5*(in[2]);
00064         }
00065 
00067         inline void 
00068         evaluateJacobian (const typename Traits::DomainType& in,         // position
00069                                           std::vector<typename Traits::JacobianType>& out) const      // return value
00070         {  
00071           out.resize(6);
00072           out[0][0][0] = sign0;       out[0][0][1] = 0;      out[0][0][2] = 0;
00073           out[0][1][0] = 0;           out[0][1][1] = 0;      out[0][1][2] = 0;
00074           out[0][2][0] = 0;           out[0][2][1] = 0;      out[0][2][2] = 0;
00075 
00076           out[1][0][0] = sign1;       out[1][0][1] = 0;      out[1][0][2] = 0; 
00077           out[1][1][0] = 0;           out[1][1][1] = 0;      out[1][1][2] = 0; 
00078           out[1][2][0] = 0;           out[1][2][1] = 0;      out[1][2][2] = 0; 
00079 
00080           out[2][0][0] = 0;           out[2][0][1] = 0;      out[2][0][2] = 0; 
00081           out[2][1][0] = 0;           out[2][1][1] = sign2;  out[2][1][2] = 0;
00082           out[2][2][0] = 0;           out[2][2][1] = 0;      out[2][2][2] = 0;
00083 
00084           out[3][0][0] = 0;           out[3][0][1] = 0;      out[3][0][2] = 0;
00085           out[3][1][0] = 0;           out[3][1][1] = sign3;  out[3][1][2] = 0;
00086           out[3][2][0] = 0;           out[3][2][1] = 0;      out[3][2][2] = 0;
00087 
00088           out[4][0][0] = 0;           out[4][0][1] = 0;      out[4][0][2] = 0;
00089           out[4][1][0] = 0;           out[4][1][1] = 0;      out[4][1][2] = 0;
00090           out[4][2][0] = 0;           out[4][2][1] = 0;      out[4][2][2] = sign4;
00091 
00092           out[5][0][0] = 0;           out[5][0][1] = 0;      out[5][0][2] = 0;
00093           out[5][1][0] = 0;           out[5][1][1] = 0;      out[5][1][2] = 0;
00094           out[5][2][0] = 0;           out[5][2][1] = 0;      out[5][2][2] = sign5;
00095         }
00096 
00098         unsigned int order () const
00099         {
00100           return 1;
00101         }
00102 
00103   private:
00104         R sign0, sign1, sign2, sign3, sign4, sign5;
00105   };
00106 
00107 
00115   template<class LB>
00116   class RT0Q3DLocalInterpolation 
00117   {
00118   public:
00119 
00121         RT0Q3DLocalInterpolation ()
00122         {
00123                 sign0 = sign1 = sign2 = sign3 = sign4 = sign5 = 1.0;
00124         }
00125 
00127         RT0Q3DLocalInterpolation (unsigned int s)
00128         {
00129           sign0 = sign1 = sign2 = sign3 = sign4 = sign5 = 1.0;
00130           if (s&1) sign0 *= -1.0;
00131           if (s&2) sign1 *= -1.0;
00132           if (s&4) sign2 *= -1.0;
00133           if (s&8) sign3 *= -1.0;
00134           if (s&16) sign4 *= -1.0;
00135           if (s&32) sign5 *= -1.0;
00136 
00137           m0[0] = 0.0; m0[1] = 0.5; m0[2] = 0.5;
00138           m1[0] = 1.0; m1[1] = 0.5; m1[2] = 0.5;
00139           m2[0] = 0.5; m2[1] = 0.0; m2[2] = 0.5;
00140           m3[0] = 0.5; m3[1] = 1.0; m3[2] = 0.5;
00141           m4[0] = 0.5; m4[1] = 0.5; m4[2] = 0.0;
00142           m5[0] = 0.5; m5[1] = 0.5; m5[2] = 1.0;
00143 
00144           n0[0] = -1.0; n0[1] =  0.0; n0[2] = 0.0;
00145           n1[0] =  1.0; n1[1] =  0.0; n1[2] = 0.0;
00146           n2[0] =  0.0; n2[1] = -1.0; n2[2] = 0.0;
00147           n3[0] =  0.0; n3[1] =  1.0; n3[2] = 0.0;
00148           n4[0] =  0.0; n4[1] =  0.0; n4[2] =-1.0;
00149           n5[0] =  0.0; n5[1] =  0.0; n5[2] = 1.0;
00150         }
00151 
00152         template<typename F, typename C>
00153         void interpolate (const F& f, std::vector<C>& out) const
00154         {
00155           // f gives v*outer normal at a point on the edge!
00156           typename F::Traits::RangeType y;
00157 
00158           out.resize(6);
00159 
00160           f.evaluate(m0,y); out[0] = (y[0]*n0[0]+y[1]*n0[1]+y[2]*n0[2])*sign0;    
00161           f.evaluate(m1,y); out[1] = (y[0]*n1[0]+y[1]*n1[1]+y[2]*n1[2])*sign1;    
00162           f.evaluate(m2,y); out[2] = (y[0]*n2[0]+y[1]*n2[1]+y[2]*n2[2])*sign2;    
00163           f.evaluate(m3,y); out[3] = (y[0]*n3[0]+y[1]*n3[1]+y[2]*n3[2])*sign3;    
00164           f.evaluate(m4,y); out[4] = (y[0]*n4[0]+y[1]*n4[1]+y[2]*n4[2])*sign4;    
00165           f.evaluate(m5,y); out[5] = (y[0]*n5[0]+y[1]*n5[1]+y[2]*n5[2])*sign5;    
00166         }
00167 
00168   private:
00169         typename LB::Traits::RangeFieldType sign0,sign1,sign2,sign3,sign4,sign5;
00170         typename LB::Traits::DomainType m0,m1,m2,m3,m4,m5;
00171         typename LB::Traits::DomainType n0,n1,n2,n3,n4,n5;
00172   };
00173 
00180   class RT0Q3DLocalCoefficients 
00181   {
00182   public:
00184         RT0Q3DLocalCoefficients () : li(6)
00185         {
00186           for (std::size_t i=0; i<6; i++)
00187                 li[i] = LocalKey(i,1,0);
00188         }
00189 
00191         std::size_t size () const
00192         {
00193           return 6;
00194         }
00195 
00197         const LocalKey& localKey (std::size_t i) const
00198         {
00199           return li[i];
00200         } 
00201 
00202   private:
00203         std::vector<LocalKey> li;
00204   };
00205 
00206 }
00207 #endif