NGSolve  4.9
fem/hcurlfe.hpp
00001 #ifndef FILE_HCURLFE
00002 #define FILE_HCURLFE
00003 
00004 /*********************************************************************/
00005 /* File:   hcurlfe.hpp                                               */
00006 /* Author: Joachim Schoeberl                                         */
00007 /* Date:   16. Apr. 2000                                             */
00008 /*********************************************************************/
00009 
00010 
00011 
00012 namespace ngfem
00013 {
00014 
00015 
00016   /*
00017     HCurl Finite Element Definitions
00018   */
00019 
00020 
00021   template <int D> class HDivFiniteElement;
00022 
00023 
00024   template <int D>
00025   class DIM_CURL_TRAIT
00026   {
00027   public:
00028     enum { DIM = (D*(D-1))/2 };
00029   };
00030 
00031 
00032   template <> class DIM_CURL_TRAIT<1>
00033   {
00034   public:
00035     enum { DIM = 0 };  // should be 0; set to 1 to make gcc34 feel ok
00036   };
00037 
00038 
00042   template <int D>
00043   class NGS_DLL_HEADER HCurlFiniteElement : public FiniteElement
00044   {
00045 
00046   public:
00047     enum { DIM = D };
00048     enum { DIM_CURL = DIM_CURL_TRAIT<D>::DIM };
00049 
00050 
00051   public:
00053     HCurlFiniteElement () { ; }
00054 
00056     HCurlFiniteElement (ELEMENT_TYPE aeltype, int andof, int aorder)
00057       : FiniteElement (aeltype, andof, aorder) { ; } 
00058   
00059     virtual ~HCurlFiniteElement () { ; }
00060 
00061     virtual string ClassName() const;
00062 
00064     virtual void CalcShape (const IntegrationPoint & ip, 
00065                             FlatMatrixFixWidth<DIM> shape) const = 0;
00066   
00068     virtual void CalcCurlShape (const IntegrationPoint & ip, 
00069                                 FlatMatrixFixWidth<DIM_CURL> curlshape) const;
00070 
00072     virtual void CalcMappedShape (const MappedIntegrationPoint<DIM,DIM> & mip,
00073                                   FlatMatrixFixWidth<DIM> shape) const;
00074 
00076     virtual void CalcMappedCurlShape (const MappedIntegrationPoint<DIM,DIM> & mip,
00077                                       FlatMatrixFixWidth<DIM_CURL> curlshape) const;
00078 
00079 
00081     const FlatMatrixFixWidth<DIM> GetShape (const IntegrationPoint & ip, 
00082                                             LocalHeap & lh) const
00083     {
00084       FlatMatrixFixWidth<DIM> shape(ndof, lh);
00085       CalcShape (ip, shape);
00086       return shape;
00087     }
00088 
00090     const FlatMatrixFixWidth<DIM_CURL> GetCurlShape (const IntegrationPoint & ip, 
00091                                                      LocalHeap & lh) const
00092     {
00093       FlatMatrixFixWidth<DIM_CURL> curlshape(ndof, lh);
00094       CalcCurlShape (ip, curlshape);
00095       return curlshape;
00096     }  
00097     
00098     template <typename TVX>
00099     Vec<DIM_CURL, typename TVX::TSCAL> 
00100     EvaluateCurlShape (const IntegrationPoint & ip, 
00101                        const TVX & x, LocalHeap & lh) const
00102     {
00103       return Trans (GetCurlShape(ip, lh)) * x;
00104     } 
00105 
00106     virtual Vec<DIM_CURL> 
00107     EvaluateCurlShape (const IntegrationPoint & ip, 
00108                        FlatVector<double> x, LocalHeap & lh) const
00109     {
00110       return Trans (GetCurlShape(ip, lh)) * x;
00111     }  
00112 
00113 
00114   protected:
00116     virtual void CalcShape1 (const IntegrationPoint & ip, 
00117                              FlatMatrixFixWidth<D> shape) const
00118     { ; }
00120     virtual void CalcShape2 (const IntegrationPoint & ip, 
00121                              FlatMatrixFixWidth<D> shape) const
00122     { ; }
00123   
00124     virtual void CalcShape3 (const IntegrationPoint & ip, 
00125                              FlatMatrixFixWidth<D> shape) const
00126     { ; }
00127   
00128     virtual void CalcShape4 (const IntegrationPoint & ip, 
00129                              FlatMatrixFixWidth<D> shape) const
00130     { ; }
00131   
00133     void ComputeEdgeMoments (int enr, ScalarFiniteElement<1> & testfe,
00134                              FlatMatrix<> moments, int order, int shape = 1) const;
00136     void ComputeFaceMoments (int fnr, HDivFiniteElement<2> & testfe,
00137                              FlatMatrix<> moments, int order, int shape = 1) const;
00139     void ComputeVolMoments (HDivFiniteElement<3> & testfe,
00140                             FlatMatrix<> moments, int order, int shape = 1) const;
00141   };
00142 
00143 
00144 
00145 
00146 
00147 
00148 
00149 
00150 
00151 
00152 
00153 
00154 
00155 
00156   // hv.DValue() = (grad u) x (grad v) 
00157   inline AutoDiff<3> Cross (const AutoDiff<3> & u,
00158                             const AutoDiff<3> & v)
00159   {
00160     AutoDiff<3> hv;
00161     hv.Value() = 0.0;
00162     hv.DValue(0) = u.DValue(1)*v.DValue(2)-u.DValue(2)*v.DValue(1);
00163     hv.DValue(1) = u.DValue(2)*v.DValue(0)-u.DValue(0)*v.DValue(2);
00164     hv.DValue(2) = u.DValue(0)*v.DValue(1)-u.DValue(1)*v.DValue(0);
00165     return hv;
00166   }
00167 
00168   inline AutoDiff<1> Cross (const AutoDiff<2> & u,
00169                             const AutoDiff<2> & v)
00170   {
00171     AutoDiff<1> hv;
00172     hv.Value() = 0.0;
00173     hv.DValue(0) = u.DValue(0)*v.DValue(1)-u.DValue(1)*v.DValue(0);
00174     return hv;
00175   }
00176 
00177 
00178 
00179 
00180 
00181 
00182 
00183   template <int DIM>
00184   class Du
00185   {
00186     enum { DIM_CURL = (DIM * (DIM-1))/2 };
00187 
00188   public:
00189     const AutoDiff<DIM> & u;
00190 
00191     Du (const AutoDiff<DIM> & au)
00192       : u(au) { ; }
00193 
00194     Vec<DIM> Value () const
00195     {
00196       Vec<DIM> val;
00197       for (int j = 0; j < DIM; j++)
00198         val(j) = u.DValue(j);
00199       return val;
00200     }
00201 
00202     Vec<DIM_CURL> CurlValue () const
00203     {
00204       return Vec<DIM> (0.0);
00205     }
00206   };
00207 
00208 
00209 
00210 
00211   template <int DIM>
00212   class uDv
00213   {
00214     enum { DIM_CURL = (DIM * (DIM-1))/2 };
00215 
00216   public:
00217     const AutoDiff<DIM> & u, v;
00218 
00219     uDv (const AutoDiff<DIM> & au, 
00220          const AutoDiff<DIM> & av)
00221       : u(au), v(av) { ; }
00222 
00223     Vec<DIM> Value () const
00224     {
00225       Vec<DIM> val;
00226       for (int j = 0; j < DIM; j++)
00227         val(j) = u.Value() * v.DValue(j);
00228       return val;
00229     }
00230 
00231     Vec<DIM_CURL> CurlValue () const
00232     {
00233       AutoDiff<DIM_CURL> hd = Cross (u, v);
00234       Vec<DIM_CURL> val;
00235       for (int i = 0; i < DIM_CURL; i++) 
00236         val(i) = hd.DValue(i);
00237       return val;
00238     }
00239   };
00240 
00241 
00242 
00243   template <int DIM>
00244   class uDv_minus_vDu
00245   {
00246     enum { DIM_CURL = (DIM * (DIM-1))/2 };
00247 
00248   public:
00249     const AutoDiff<DIM> & u, v;
00250 
00251     uDv_minus_vDu (const AutoDiff<DIM> & au, 
00252                    const AutoDiff<DIM> & av)
00253       : u(au), v(av) { ; }
00254 
00255     Vec<DIM> Value () const
00256     {
00257       Vec<DIM> val;
00258       for (int j = 0; j < DIM; j++)
00259         val(j) = u.Value() * v.DValue(j) - v.Value() * u.DValue(j);
00260       return val;
00261     }
00262 
00263     Vec<DIM_CURL> CurlValue () const
00264     {
00265       AutoDiff<DIM_CURL> hd = Cross (u, v);
00266       Vec<DIM_CURL> val;
00267       for (int i = 0; i < DIM_CURL; i++) 
00268         val(i) = 2 * hd.DValue(i);
00269       return val;
00270     }
00271   };
00272 
00273 
00274   template <int DIM>
00275   class wuDv_minus_wvDu
00276   {
00277     enum { DIM_CURL = (DIM * (DIM-1))/2 };
00278 
00279   public:
00280     const AutoDiff<DIM> & u, v, w;
00281 
00282     wuDv_minus_wvDu (const AutoDiff<DIM> & au, 
00283                      const AutoDiff<DIM> & av,
00284                      const AutoDiff<DIM> & aw)
00285       : u(au), v(av), w(aw) { ; }
00286 
00287     Vec<DIM> Value () const
00288     {
00289       Vec<DIM> val;
00290       for (int j = 0; j < DIM; j++)
00291         val(j) = w.Value() * (u.Value() * v.DValue(j) - v.Value() * u.DValue(j));
00292       return val;
00293     }
00294 
00295     Vec<DIM_CURL> CurlValue () const
00296     {
00297       AutoDiff<DIM_CURL> hd = Cross (u*w, v) + Cross(u, v*w);
00298       Vec<DIM_CURL> val;
00299       for (int i = 0; i < DIM_CURL; i++) 
00300         val(i) = hd.DValue(i);
00301       return val;
00302     }
00303   };
00304 
00305 
00306 
00307 
00308 
00309 
00310   template <int DIM>
00311   class HCurlShapeElement
00312   {
00313     FlatVec<DIM> data;
00314   public:
00315     HCurlShapeElement (FlatVec<DIM> adata) : data(adata) { ; }
00316 
00317     template <typename F1>
00318     void operator= (F1 form1)  { data = form1.Value(); }
00319   };
00320 
00321   template <int DIM>
00322   class HCurlCurlShapeElement
00323   {
00324     enum { DIM_CURL = (DIM * (DIM-1))/2 };
00325     FlatVec<DIM_CURL> data;
00326   public:
00327     HCurlCurlShapeElement (FlatVec<DIM_CURL> adata) : data(adata) { ; }
00328 
00329     template <typename F1>
00330     void operator= (F1 form1)  { data = form1.CurlValue(); }
00331   };
00332 
00333 
00334   template <int DIM>
00335   class HCurlEvaluateCurlElement
00336   {
00337     const double & coef;
00338     enum { DIM_CURL = (DIM * (DIM-1))/2 };
00339     Vec<DIM_CURL> & sum;
00340   public:
00341     HCurlEvaluateCurlElement (const double & acoef, Vec<DIM_CURL> & asum)
00342       : coef(acoef), sum(asum) { ; }
00343 
00344     template <typename F1>
00345     void operator= (F1 form1)  { sum += coef * form1.CurlValue(); }
00346   };
00347 
00348 
00349 
00350   template <int DIM>
00351   class HCurlShapeAssign
00352   {
00353     FlatMatrixFixWidth<DIM> shape; 
00354   public:
00355     HCurlShapeAssign (FlatMatrixFixWidth<DIM> mat) : shape(mat) { ; }
00356 
00357     HCurlShapeElement<DIM> operator[] (int i) const
00358     { return HCurlShapeElement<DIM> (shape.Row(i)); }
00359   };
00360 
00361   template <int DIM>
00362   class HCurlCurlShapeAssign
00363   {
00364     enum { DIM_CURL = (DIM * (DIM-1))/2 };
00365     FlatMatrixFixWidth<DIM_CURL> cshape; 
00366   public:
00367     HCurlCurlShapeAssign (FlatMatrixFixWidth<DIM_CURL> mat) : cshape(mat) { ; }
00368 
00369     HCurlCurlShapeElement<DIM> operator[] (int i) const 
00370     { return HCurlCurlShapeElement<DIM> (cshape.Row(i)); }
00371   };
00372 
00373   template <int DIM>
00374   class HCurlEvaluateCurl
00375   {
00376     FlatVector<> coefs;
00377     enum { DIM_CURL = (DIM * (DIM-1))/2 };
00378     Vec<DIM_CURL> sum;
00379   public:
00380     HCurlEvaluateCurl (FlatVector<> acoefs) : coefs(acoefs), sum(0.0) { ; }
00381 
00382     HCurlEvaluateCurlElement<DIM> operator[] (int i) 
00383     { return HCurlEvaluateCurlElement<DIM> (coefs(i), sum); }
00384 
00385     Vec<DIM_CURL> Sum() { return sum; }
00386   };
00387 
00388 
00389 
00390 
00391 
00392 
00393 
00399   template <class FEL, ELEMENT_TYPE ET, int NDOF, int ORDER>
00400   class T_HCurlFiniteElement : public HCurlFiniteElement<ET_trait<ET>::DIM>
00401   {
00402 
00403   public:
00404     
00405   protected:
00406     enum { DIM = ET_trait<ET>::DIM };
00407     using HCurlFiniteElement<DIM>::DIM_CURL;
00408 
00409     T_HCurlFiniteElement ()
00410       : HCurlFiniteElement<DIM> (ET, NDOF, ORDER) { ; }
00411 
00412     virtual ~T_HCurlFiniteElement() { ; }
00413 
00414   public:
00415 
00416     virtual void CalcShape (const IntegrationPoint & ip, 
00417                             FlatMatrixFixWidth<DIM> shape) const
00418     {
00419       AutoDiff<DIM> adp[DIM];
00420       for (int i = 0; i < DIM; i++)
00421         adp[i] = AutoDiff<DIM> (ip(i), i);
00422       
00423       HCurlShapeAssign<DIM> ds(shape); 
00424       FEL::T_CalcShape (adp, ds);
00425     }
00426 
00427     virtual void
00428     CalcMappedShape (const MappedIntegrationPoint<DIM,DIM> & mip,
00429                      FlatMatrixFixWidth<DIM> shape) const
00430     {
00431       AutoDiff<DIM> adp[DIM];
00432       
00433       for (int i = 0; i < DIM; i++)
00434         adp[i].Value() = mip.IP()(i);
00435       
00436       for (int i = 0; i < DIM; i++)
00437         for (int j = 0; j < DIM; j++)
00438           adp[i].DValue(j) = mip.GetJacobianInverse()(i,j);
00439       
00440       HCurlShapeAssign<DIM> ds(shape); 
00441       FEL::T_CalcShape (adp, ds);
00442     }
00443 
00444 
00445     virtual void CalcCurlShape (const IntegrationPoint & ip, 
00446                                 FlatMatrixFixWidth<DIM_CURL> curlshape) const
00447     {
00448       AutoDiff<DIM> adp[DIM];
00449       for (int i = 0; i < DIM; i++)
00450         adp[i] = AutoDiff<DIM> (ip(i), i);
00451 
00452       HCurlCurlShapeAssign<DIM> ds(curlshape); 
00453       FEL::T_CalcShape (adp, ds);
00454     }
00455 
00456     virtual void
00457     CalcMappedCurlShape (const MappedIntegrationPoint<DIM,DIM> & mip,
00458                          FlatMatrixFixWidth<DIM_CURL> curlshape) const
00459     {
00460       AutoDiff<DIM> adp[DIM];
00461 
00462       for (int i = 0; i < DIM; i++)
00463         adp[i].Value() = mip.IP()(i);
00464       
00465       for (int i = 0; i < DIM; i++)
00466         for (int j = 0; j < DIM; j++)
00467           adp[i].DValue(j) = mip.GetJacobianInverse()(i,j);
00468 
00469       HCurlCurlShapeAssign<DIM> ds(curlshape); 
00470       FEL::T_CalcShape (adp, ds);
00471     }
00472 
00473     virtual Vec <DIM_CURL>
00474     EvaluateCurlShape (const IntegrationPoint & ip, 
00475                        FlatVector<double> x,
00476                        LocalHeap & lh) const
00477     {
00478       AutoDiff<DIM> adp[DIM];
00479       for (int i = 0; i < DIM; i++)
00480         adp[i] = AutoDiff<DIM> (ip(i), i);
00481       
00482       HCurlEvaluateCurl<DIM> ds(x); 
00483       FEL::T_CalcShape (adp, ds);
00484       return ds.Sum();
00485     }
00486 
00487   };
00488 
00489 
00490 
00491 
00492   template <int D>
00493   extern void ComputeGradientMatrix (const ScalarFiniteElement<D> & h1fe,
00494                                      const HCurlFiniteElement<D> & hcurlfe,
00495                                      FlatMatrix<> gradient);
00496 
00497 }
00498 
00499 
00500 
00501 
00502 #endif