IPPL: src/expde/formulas/diffopv.h Source File

00001 //    expde: expression templates for partial differential equations.
00002 //    Copyright (C) 2001  Christoph Pflaum
00003 //    This program is free software; you can redistribute it and/or modify
00004 //    it under the terms of the GNU General Public License as published by
00005 //    the Free Software Foundation; either version 2 of the License, or
00006 //    (at your option) any later version.
00007 //
00008 //    This program is distributed in the hope that it will be useful,
00009 //    but WITHOUT ANY WARRANTY; without even the implied warranty of
00010 //    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00011 //    GNU General Public License for more details.
00012 //
00013 //                 SEE  Notice1.doc made by 
00014 //                 LAWRENCE LIVERMORE NATIONAL LABORATORY
00015 //
00016 
00017 // ------------------------------------------------------------
00018 // diffopv.h
00019 //
00020 // ------------------------------------------------------------
00021 
00022 #ifndef DIFFOPV_H_
00023 #define DIFFOPV_H_
00024         
00025 // -----------------------------------------------------------
00026 //    FFFF    OO     RR     M     M   EEEE  L     N     N
00027 //    F      O  O    R R    M M M M   E     L     N N   N
00028 //    FFF    O  O    RR     M  M  M   EEEE  L     N  N  N
00029 //    F      O  O    R R    M     M   E     L     N   N N
00030 //    F       OO     R  R   M     M   EEEE  LLLL  N     N
00031 // -----------------------------------------------------------
00032 
00034 // 4.:  Variable coefficient Differential operators
00036 
00037 #define var_interpo_EST (aEST+aM)
00038 #define var_interpo_WST (aST+aW)
00039 #define var_interpo_ENT (aET+aN)
00040 #define var_interpo_ESD (aES+aD)
00041 #define var_interpo_WND (aWND+aM)
00042 #define var_interpo_END (aND+aE)
00043 #define var_interpo_WSD (aWD+aS)
00044 #define var_interpo_WNT (aWN+aT)
00045 
00046 
00047 #define Coefficient_triangle  (bocedata->vars[tet->N0()][num_coefficient]+bocedata->vars[tet->N1()][num_coefficient]+bocedata->vars[tet->N2()][num_coefficient]+bocedata->vars[tet->N3()][num_coefficient])/4.0
00048 
00049 
00050 // Formeln fuer den Laplace Operator
00051 
00052 class laplace_FE_variable {
00053  public:
00054   diff_type typ_u() { return  Poisson_type; }; 
00055   diff_type typ_v() { return  Poisson_type; };
00056 
00057   laplace_FE_variable() { }
00058   static inline int    ops_interior() { return 68; }
00059   static inline int    ops_diag_interior() { return 20; }
00060   static inline double stencil(double uN,
00061                                double uW, double uM, double uE,
00062                                double uS, double uT, double uD,
00063                                double uND, double uWN, double uWT,
00064                                double uED, double uST, double uES,
00065                                double uEST, double uWND,
00066                                double aN,
00067                                double aW, double aM, double aE,
00068                                double aS, double aT, double aD,
00069                                double aST, double aET, double aES,
00070                                double aND, double aWD, double aWN,
00071                                double aEST, double aWND,
00072                                double h_mesh) {
00073     return  (var_interpo_WSD * (3.0*uM - uW - uS - uD) +
00074              var_interpo_ESD * (5.0*uM - 3.0*uD - uE - uS - uES + uED) +
00075              var_interpo_WND * (7.0*uM - 3.0*(uW+uD) -uN-uND-uWN + 2.0*uWND) +
00076              var_interpo_END * (5.0*uM - 3.0*uE - uN - uD - uND + uED) +
00077              var_interpo_WST * (5.0*uM - 3.0*uW - uT - uS - uST + uWT) +
00078              var_interpo_EST * (7.0*uM - 3.0*(uE+uT) -uS-uST-uES + 2.0*uEST) +
00079              var_interpo_WNT * (5.0*uM - 3.0*uT - uN - uW - uWN + uWT) +
00080              var_interpo_ENT * (3.0*uM - uE - uN - uT)
00081              ) * h_mesh / 12.0; // additional divide by 2.0
00082   }
00083   static inline double diag_stencil(double h_mesh,
00084                                     double aN,
00085                                     double aW, double aM, double aE,
00086                                     double aS, double aT, double aD,
00087                                     double aST, double aET, double aES,
00088                                     double aND, double aWD, double aWN,
00089                                     double aEST, double aWND) {
00090     return ((var_interpo_WSD + var_interpo_ENT) * 3.0 +
00091             (var_interpo_ESD + var_interpo_END +
00092              var_interpo_WST + var_interpo_WNT) * 5.0 +
00093             (var_interpo_WND + var_interpo_EST) * 7.0) / 12.0 * h_mesh;
00094   }
00095   static inline double Factor_reg_cell(double h_mesh) {
00096     return h_mesh / 6.0;
00097   }
00098   static inline double diag_F_reg_cell(dir_sons dir_v,
00099                                        double h_mesh, double* sten) {
00100     return sten[dir_v*9] * h_mesh / 6.0;
00101   }
00102 
00103   // Operator applied to variable with number: num_variable
00104   // with test function: num_v
00105   // on a boundary cell
00106   static inline double F_bo_cell(const BoCeData* bocedata,
00107                                  int num_v, int num_variable,
00108                                  int num_coefficient) {
00109     static double sum;
00110     Tetraeder_storage* tet;
00111     
00112     sum=0.0;
00113     for(tet=bocedata->Give_tets();tet!=NULL;tet=tet->Next()) {
00114       if(tet->N0()==num_v) {
00115         sum = sum + (
00116                      XM0*(XM0*UU0 + XM1*UU1 + XM2*UU2 + XM3*UU3) +
00117                      YM0*(YM0*UU0 + YM1*UU1 + YM2*UU2 + YM3*UU3) +
00118                      ZM0*(ZM0*UU0 + ZM1*UU1 + ZM2*UU2 + ZM3*UU3)
00119                      ) / (tet->Det()*6.0)
00120           * Coefficient_triangle;
00121       }
00122       else if(tet->N1()==num_v) {
00123         sum = sum + (
00124                      XM1*(XM0*UU0 + XM1*UU1 + XM2*UU2 + XM3*UU3) +
00125                      YM1*(YM0*UU0 + YM1*UU1 + YM2*UU2 + YM3*UU3) +
00126                      ZM1*(ZM0*UU0 + ZM1*UU1 + ZM2*UU2 + ZM3*UU3)
00127                      ) / (tet->Det()*6.0)
00128           * Coefficient_triangle;
00129       }
00130       else if(tet->N2()==num_v) {
00131         sum = sum + (
00132                      XM2*(XM0*UU0 + XM1*UU1 + XM2*UU2 + XM3*UU3) +
00133                      YM2*(YM0*UU0 + YM1*UU1 + YM2*UU2 + YM3*UU3) +
00134                      ZM2*(ZM0*UU0 + ZM1*UU1 + ZM2*UU2 + ZM3*UU3)
00135                      ) / (tet->Det()*6.0)
00136           * Coefficient_triangle;
00137       }
00138       else if(tet->N3()==num_v) {
00139         sum = sum + (
00140                      XM3*(XM0*UU0 + XM1*UU1 + XM2*UU2 + XM3*UU3) +
00141                      YM3*(YM0*UU0 + YM1*UU1 + YM2*UU2 + YM3*UU3) +
00142                      ZM3*(ZM0*UU0 + ZM1*UU1 + ZM2*UU2 + ZM3*UU3)
00143                      ) / (tet->Det()*6.0)
00144           * Coefficient_triangle;
00145       }
00146     }
00147     return sum;
00148   }
00149 
00150   // Operator applied to a vector U_array
00151   // with test function: num_v
00152   // on a boundary cell
00153   static inline double Local_matrix_bo_cell(const BoCeData* bocedata,
00154                                             double*  U_array, int num_v,
00155                                             int num_coefficient) {
00156     static double sum;
00157     Tetraeder_storage* tet;
00158     
00159     sum=0.0;
00160     for(tet=bocedata->Give_tets();tet!=NULL;tet=tet->Next()) {
00161       if(tet->N0()==num_v) {
00162         sum = sum + (
00163                      XM0*(XM0*UAA0 + XM1*UAA1 + XM2*UAA2 + XM3*UAA3) +
00164                      YM0*(YM0*UAA0 + YM1*UAA1 + YM2*UAA2 + YM3*UAA3) +
00165                      ZM0*(ZM0*UAA0 + ZM1*UAA1 + ZM2*UAA2 + ZM3*UAA3)) /
00166           (tet->Det()*6.0) * Coefficient_triangle;
00167       }
00168       else if(tet->N1()==num_v) {
00169         sum = sum + (
00170                      XM1*(XM0*UAA0 + XM1*UAA1 + XM2*UAA2 + XM3*UAA3) +
00171                      YM1*(YM0*UAA0 + YM1*UAA1 + YM2*UAA2 + YM3*UAA3) +
00172                      ZM1*(ZM0*UAA0 + ZM1*UAA1 + ZM2*UAA2 + ZM3*UAA3)) /
00173           (tet->Det()*6.0) * Coefficient_triangle;
00174       }
00175       else if(tet->N2()==num_v) {
00176         sum = sum + (
00177                      XM2*(XM0*UAA0 + XM1*UAA1 + XM2*UAA2 + XM3*UAA3) +
00178                      YM2*(YM0*UAA0 + YM1*UAA1 + YM2*UAA2 + YM3*UAA3) +
00179                      ZM2*(ZM0*UAA0 + ZM1*UAA1 + ZM2*UAA2 + ZM3*UAA3)) /
00180           (tet->Det()*6.0) * Coefficient_triangle;
00181       }
00182       else if(tet->N3()==num_v) {
00183         sum = sum + (
00184                      XM3*(XM0*UAA0 + XM1*UAA1 + XM2*UAA2 + XM3*UAA3) +
00185                      YM3*(YM0*UAA0 + YM1*UAA1 + YM2*UAA2 + YM3*UAA3) +
00186                      ZM3*(ZM0*UAA0 + ZM1*UAA1 + ZM2*UAA2 + ZM3*UAA3)) /
00187           (tet->Det()*6.0) * Coefficient_triangle;
00188       }
00189     }
00190     return sum;  
00191   }
00192 
00193   static inline double diag_F_bo_cell(const BoCeData* bocedata,
00194                                       int num_v, int num_coefficient) {
00195     static double sum;
00196     Tetraeder_storage* tet;
00197     
00198     sum=0.0;
00199     
00200     for(tet=bocedata->Give_tets();tet!=NULL;tet=tet->Next()) {
00201       if(tet->N0()==num_v) {
00202         sum = sum + (XM0*XM0 + YM0*YM0 + ZM0*ZM0) / (tet->Det()*6.0)
00203           * Coefficient_triangle;
00204       }
00205       else if(tet->N1()==num_v) {
00206         sum = sum + (XM1*XM1 + YM1*YM1 + ZM1*ZM1) / (tet->Det()*6.0)
00207           * Coefficient_triangle;
00208       }
00209       else if(tet->N2()==num_v) {
00210         sum = sum + (XM2*XM2 + YM2*YM2 + ZM2*ZM2) / (tet->Det()*6.0)
00211           * Coefficient_triangle;
00212       }
00213       else if(tet->N3()==num_v) {
00214         sum = sum + (XM3*XM3 + YM3*YM3 + ZM3*ZM3) / (tet->Det()*6.0)
00215           * Coefficient_triangle;
00216       }
00217     }
00218     return sum;
00219   }
00220 };
00221 
00222 
00223 // Formeln fuer dxdx
00224 class dxdx_FE_variable {
00225  public:
00226   diff_type typ_u() { return  Dx_type; }; 
00227   diff_type typ_v() { return  Dx_type; };
00228 
00229   dxdx_FE_variable() { }
00230   static inline int    ops_interior() { return 34; }
00231   static inline int    ops_diag_interior() { return 18; }
00232   static inline double stencil(double uN,
00233                                double uW, double uM, double uE,
00234                                double uS, double uT, double uD,
00235                                double uND, double uWN, double uWT,
00236                                double uED, double uST, double uES,
00237                                double uEST, double uWND,
00238                                double aN,
00239                                double aW, double aM, double aE,
00240                                double aS, double aT, double aD,
00241                                double aST, double aET, double aES,
00242                                double aND, double aWD, double aWN,
00243                                double aEST, double aWND,
00244                                double h_mesh) {
00245     return  (var_interpo_WSD * (uM - uW) +
00246              var_interpo_ESD * (uM - uE) +
00247              (var_interpo_WND * (uM - uW) +
00248               var_interpo_END * (uM - uE) +
00249               var_interpo_WST * (uM - uW) +
00250               var_interpo_EST * (uM - uE))* 2.0 +
00251              var_interpo_WNT * (uM - uW) +
00252              var_interpo_ENT * (uM - uE) 
00253              ) * h_mesh / 12.0; // additional divide by 2.0
00254   }
00255   static inline double diag_stencil(double h_mesh,
00256                                     double aN,
00257                                     double aW, double aM, double aE,
00258                                     double aS, double aT, double aD,
00259                                     double aST, double aET, double aES,
00260                                     double aND, double aWD, double aWN,
00261                                     double aEST, double aWND) {
00262     return  (var_interpo_WSD +
00263              var_interpo_ESD +
00264              (var_interpo_WND +
00265               var_interpo_END +
00266               var_interpo_WST +
00267               var_interpo_EST) * 2.0 +
00268              var_interpo_WNT +
00269              var_interpo_ENT 
00270              ) * h_mesh / 12.0;
00271   }
00272   static inline double Factor_reg_cell(double h_mesh) {
00273     return h_mesh / 6.0;
00274   }
00275   static inline double diag_F_reg_cell(dir_sons dir_v,
00276                                        double h_mesh, double* sten) {
00277     return sten[dir_v*9] * h_mesh / 6.0;
00278   }
00279 
00280   // Operator applied to variable with number
00281   static inline double F_bo_cell(const BoCeData* bocedata,
00282                                  int num_v, int num_variable,
00283                                  int num_coefficient) {
00284     static double sum;
00285     Tetraeder_storage* tet;
00286 
00287     sum=0.0;
00288     for(tet=bocedata->Give_tets();tet!=NULL;tet=tet->Next()) {
00289       if(tet->N0()==num_v) {
00290         sum = sum +
00291               XM0*(XM0*UU0 + XM1*UU1 + XM2*UU2 + XM3*UU3) / (tet->Det()*6.0)
00292           * Coefficient_triangle;
00293       }
00294       else if(tet->N1()==num_v) {
00295         sum = sum +
00296               XM1*(XM0*UU0 + XM1*UU1 + XM2*UU2 + XM3*UU3) / (tet->Det()*6.0)
00297           * Coefficient_triangle;
00298       }
00299       else if(tet->N2()==num_v) {
00300         sum = sum +
00301               XM2*(XM0*UU0 + XM1*UU1 + XM2*UU2 + XM3*UU3) / (tet->Det()*6.0)
00302           * Coefficient_triangle;
00303       }
00304       else if(tet->N3()==num_v) {
00305         sum = sum +
00306               XM3*(XM0*UU0 + XM1*UU1 + XM2*UU2 + XM3*UU3) / (tet->Det()*6.0)
00307           * Coefficient_triangle;
00308       }
00309     }    
00310     return sum;
00311   }
00312 
00313   // Operator applied to a vector U_array
00314   static inline double Local_matrix_bo_cell(const BoCeData* bocedata,
00315                                             double*  U_array, int num_v,
00316                                             int num_coefficient) {
00317     static double sum;
00318     Tetraeder_storage* tet;
00319 
00320     sum=0.0;
00321     for(tet=bocedata->Give_tets();tet!=NULL;tet=tet->Next()) {
00322       if(tet->N0()==num_v) {
00323         sum = sum +
00324           XM0*(XM0*UAA0 + XM1*UAA1 + XM2*UAA2 + XM3*UAA3) / (tet->Det()*6.0)
00325           * Coefficient_triangle;
00326       }
00327       else if(tet->N1()==num_v) {
00328         sum = sum +
00329           XM1*(XM0*UAA0 + XM1*UAA1 + XM2*UAA2 + XM3*UAA3) / (tet->Det()*6.0)
00330           * Coefficient_triangle;
00331       }
00332       else if(tet->N2()==num_v) {
00333         sum = sum +
00334           XM2*(XM0*UAA0 + XM1*UAA1 + XM2*UAA2 + XM3*UAA3) / (tet->Det()*6.0)
00335           * Coefficient_triangle;
00336       }
00337       else if(tet->N3()==num_v) {
00338         sum = sum +
00339           XM3*(XM0*UAA0 + XM1*UAA1 + XM2*UAA2 + XM3*UAA3) / (tet->Det()*6.0)
00340           * Coefficient_triangle;
00341       }
00342     }    
00343     return sum;
00344   }
00345 
00346   static inline double diag_F_bo_cell(const BoCeData* bocedata,
00347                                       int num_v, int num_coefficient) {
00348    static double sum;
00349     Tetraeder_storage* tet;
00350 
00351     sum=0.0;
00352     for(tet=bocedata->Give_tets();tet!=NULL;tet=tet->Next()) {
00353       if(tet->N0()==num_v) {
00354         sum = sum + XM0*XM0 / (tet->Det()*6.0)
00355           * Coefficient_triangle;
00356       }
00357       else if(tet->N1()==num_v) {
00358         sum = sum + XM1*XM1 / (tet->Det()*6.0)
00359           * Coefficient_triangle;
00360       }
00361       else if(tet->N2()==num_v) {
00362         sum = sum + XM2*XM2 / (tet->Det()*6.0)
00363           * Coefficient_triangle;
00364       }
00365       else if(tet->N3()==num_v) {
00366         sum = sum + XM3*XM3 / (tet->Det()*6.0)
00367           * Coefficient_triangle;
00368       }
00369     }     
00370     return sum;
00371   }
00372 };
00373 
00374 
00375 
00376 // Formeln fuer dydy
00377 class dydy_FE_variable {
00378  public:
00379   diff_type typ_u() { return  Dy_type; }; 
00380   diff_type typ_v() { return  Dy_type; };
00381 
00382   dydy_FE_variable() { }
00383   static inline int    ops_interior() { return 63; } // 55+8
00384   static inline int    ops_diag_interior() { return 19; }
00385   static inline double stencil(double uN,
00386                                double uW, double uM, double uE,
00387                                double uS, double uT, double uD,
00388                                double uND, double uWN, double uWT,
00389                                double uED, double uST, double uES,
00390                                double uEST, double uWND,
00391                                double aN,
00392                                double aW, double aM, double aE,
00393                                double aS, double aT, double aD,
00394                                double aST, double aET, double aES,
00395                                double aND, double aWD, double aWN,
00396                                double aEST, double aWND,
00397                                double h_mesh) {
00398     return  (var_interpo_WSD * (uM - uS) +
00399              var_interpo_ESD * (2.0*uM - uS - uES - uD + uED) +
00400              var_interpo_WND * (3.0*uM - uN - uW - uD  - uWN - uND + 2.0*uWND)+
00401              var_interpo_END * (2.0*uM - uE - uN - uND + uED) +
00402              var_interpo_WST * (2.0*uM - uW - uS - uST + uWT) +
00403              var_interpo_EST * (3.0*uM - uT - uE - uS  - uST - uES + 2.0*uEST)+
00404              var_interpo_WNT * (2.0*uM - uT - uN - uWN + uWT) +
00405              var_interpo_ENT * (uM - uN) 
00406              ) * h_mesh / 12.0; // additional divide by 2.0
00407   }
00408   static inline double diag_stencil(double h_mesh,
00409                                     double aN,
00410                                     double aW, double aM, double aE,
00411                                     double aS, double aT, double aD,
00412                                     double aST, double aET, double aES,
00413                                     double aND, double aWD, double aWN,
00414                                     double aEST, double aWND) {
00415     return  (var_interpo_WSD +
00416              var_interpo_ENT +
00417              (var_interpo_WND +
00418               var_interpo_EST) * 3.0 +
00419              (var_interpo_END +
00420               var_interpo_WST +
00421               var_interpo_ESD +
00422               var_interpo_WNT) * 2.0
00423              ) * h_mesh / 12.0; // additional divide by 2.0
00424   }
00425   static inline double Factor_reg_cell(double h_mesh) {
00426     return h_mesh / 6.0;
00427   }
00428   static inline double diag_F_reg_cell(dir_sons dir_v,
00429                                        double h_mesh, double* sten) {
00430     return sten[dir_v*9] * h_mesh / 6.0;
00431   }
00432 
00433   // Operator applied to variable with number
00434   static inline double F_bo_cell(const BoCeData* bocedata,
00435                                  int num_v, int num_variable,
00436                                  int num_coefficient) {
00437     static double sum;
00438     Tetraeder_storage* tet;
00439 
00440     sum=0.0;
00441     for(tet=bocedata->Give_tets();tet!=NULL;tet=tet->Next()) {
00442       if(tet->N0()==num_v) {
00443         sum = sum +
00444               YM0*(YM0*UU0 + YM1*UU1 + YM2*UU2 + YM3*UU3) / (tet->Det()*6.0)
00445           * Coefficient_triangle;
00446       }
00447       else if(tet->N1()==num_v) {
00448         sum = sum +
00449               YM1*(YM0*UU0 + YM1*UU1 + YM2*UU2 + YM3*UU3) / (tet->Det()*6.0)
00450           * Coefficient_triangle;
00451       }
00452       else if(tet->N2()==num_v) {
00453         sum = sum +
00454               YM2*(YM0*UU0 + YM1*UU1 + YM2*UU2 + YM3*UU3) / (tet->Det()*6.0)
00455           * Coefficient_triangle;
00456       }
00457       else if(tet->N3()==num_v) {
00458         sum = sum +
00459               YM3*(YM0*UU0 + YM1*UU1 + YM2*UU2 + YM3*UU3) / (tet->Det()*6.0)
00460           * Coefficient_triangle;
00461       }
00462     }    
00463     return sum;
00464   }
00465   // Operator applied to a vector U_array
00466   static inline double Local_matrix_bo_cell(const BoCeData* bocedata,
00467                                             double*  U_array, int num_v,
00468                                             int num_coefficient) {
00469     static double sum;
00470     Tetraeder_storage* tet;
00471 
00472     sum=0.0;
00473     for(tet=bocedata->Give_tets();tet!=NULL;tet=tet->Next()) {
00474       if(tet->N0()==num_v) {
00475         sum = sum +
00476           YM0*(YM0*UAA0 + YM1*UAA1 + YM2*UAA2 + YM3*UAA3) / (tet->Det()*6.0)
00477           * Coefficient_triangle;
00478       }
00479       else if(tet->N1()==num_v) {
00480         sum = sum +
00481           YM1*(YM0*UAA0 + YM1*UAA1 + YM2*UAA2 + YM3*UAA3) / (tet->Det()*6.0)
00482           * Coefficient_triangle;
00483       }
00484       else if(tet->N2()==num_v) {
00485         sum = sum +
00486           YM2*(YM0*UAA0 + YM1*UAA1 + YM2*UAA2 + YM3*UAA3) / (tet->Det()*6.0)
00487           * Coefficient_triangle;
00488       }
00489       else if(tet->N3()==num_v) {
00490         sum = sum +
00491           YM3*(YM0*UAA0 + YM1*UAA1 + YM2*UAA2 + YM3*UAA3) / (tet->Det()*6.0)
00492           * Coefficient_triangle;
00493       }
00494     }    
00495     return sum;
00496   }
00497   static inline double diag_F_bo_cell(const BoCeData* bocedata,
00498                                       int num_v, int num_coefficient) {
00499    static double sum;
00500     Tetraeder_storage* tet;
00501 
00502     sum=0.0;
00503     for(tet=bocedata->Give_tets();tet!=NULL;tet=tet->Next()) {
00504       if(tet->N0()==num_v) {
00505         sum = sum + YM0*YM0 / (tet->Det()*6.0)
00506           * Coefficient_triangle;
00507       }
00508       else if(tet->N1()==num_v) {
00509         sum = sum + YM1*YM1 / (tet->Det()*6.0)
00510           * Coefficient_triangle;
00511       }
00512       else if(tet->N2()==num_v) {
00513         sum = sum + YM2*YM2 / (tet->Det()*6.0)
00514           * Coefficient_triangle;
00515       }
00516       else if(tet->N3()==num_v) {
00517         sum = sum + YM3*YM3 / (tet->Det()*6.0)
00518           * Coefficient_triangle;
00519       }
00520     }     
00521     return sum;
00522   }
00523 };
00524 
00525 
00526 // Formeln fuer dzdz
00527 class dzdz_FE_variable {
00528  public:  
00529   diff_type typ_u() { return  Dz_type; }; 
00530   diff_type typ_v() { return  Dz_type; };
00531 
00532   dzdz_FE_variable() { }
00533   static inline int    ops_interior() { return 34; } // 26 + 8
00534   static inline int    ops_diag_interior() { return 18; }
00535   static inline double stencil(double uN,
00536                                double uW, double uM, double uE,
00537                                double uS, double uT, double uD,
00538                                double uND, double uWN, double uWT,
00539                                double uED, double uST, double uES,
00540                                double uEST, double uWND,
00541                                double aN,
00542                                double aW, double aM, double aE,
00543                                double aS, double aT, double aD,
00544                                double aST, double aET, double aES,
00545                                double aND, double aWD, double aWN,
00546                                double aEST, double aWND,
00547                                double h_mesh) {
00548     return  (var_interpo_WSD * (uM - uD) +
00549              var_interpo_END * (uM - uD) +
00550              var_interpo_WST * (uM - uT) +
00551              (var_interpo_WND * (uM - uD) +
00552               var_interpo_ESD * (uM - uD) +
00553               var_interpo_WNT * (uM - uT) +
00554               var_interpo_EST * (uM - uT)) * 2.0 +
00555              var_interpo_ENT * (uM - uT) 
00556              ) * h_mesh / 12.0; // additional divide by 2.0
00557   }
00558   static inline double diag_stencil(double h_mesh,
00559                                     double aN,
00560                                     double aW, double aM, double aE,
00561                                     double aS, double aT, double aD,
00562                                     double aST, double aET, double aES,
00563                                     double aND, double aWD, double aWN,
00564                                     double aEST, double aWND) {
00565     return  (var_interpo_WSD +
00566              var_interpo_END +
00567              var_interpo_WST +
00568              (var_interpo_WND +
00569               var_interpo_ESD +
00570               var_interpo_WNT +
00571               var_interpo_EST) * 2.0 +
00572              var_interpo_ENT 
00573              ) * h_mesh / 12.0; // additional divide by 2.0
00574   }
00575   static inline double Factor_reg_cell(double h_mesh) {
00576     return h_mesh / 6.0;
00577   }
00578   static inline double diag_F_reg_cell(dir_sons dir_v,
00579                                        double h_mesh, double* sten) {
00580     return sten[dir_v*9] * h_mesh / 6.0;
00581   }
00582 
00583   // Operator applied to variable with number
00584   static inline double F_bo_cell(const BoCeData* bocedata,
00585                                  int num_v, int num_variable,
00586                                  int num_coefficient) {
00587     static double sum;
00588     Tetraeder_storage* tet;
00589 
00590     sum=0.0;
00591     for(tet=bocedata->Give_tets();tet!=NULL;tet=tet->Next()) {
00592       if(tet->N0()==num_v) {
00593         sum = sum +
00594               ZM0*(ZM0*UU0 + ZM1*UU1 + ZM2*UU2 + ZM3*UU3) / (tet->Det()*6.0)
00595           * Coefficient_triangle;
00596       }
00597       else if(tet->N1()==num_v) {
00598         sum = sum +
00599               ZM1*(ZM0*UU0 + ZM1*UU1 + ZM2*UU2 + ZM3*UU3) / (tet->Det()*6.0)
00600           * Coefficient_triangle;
00601       }
00602       else if(tet->N2()==num_v) {
00603         sum = sum +
00604               ZM2*(ZM0*UU0 + ZM1*UU1 + ZM2*UU2 + ZM3*UU3) / (tet->Det()*6.0)
00605           * Coefficient_triangle;
00606       }
00607       else if(tet->N3()==num_v) {
00608         sum = sum +
00609               ZM3*(ZM0*UU0 + ZM1*UU1 + ZM2*UU2 + ZM3*UU3) / (tet->Det()*6.0)
00610           * Coefficient_triangle;
00611       }
00612     }     
00613     return sum;
00614   }
00615   // Operator applied to a vector U_array
00616   static inline double Local_matrix_bo_cell(const BoCeData* bocedata,
00617                                             double*  U_array, int num_v,
00618                                             int num_coefficient) {
00619     static double sum;
00620     Tetraeder_storage* tet;
00621 
00622     sum=0.0;
00623     for(tet=bocedata->Give_tets();tet!=NULL;tet=tet->Next()) {
00624       if(tet->N0()==num_v) {
00625         sum = sum +
00626           ZM0*(ZM0*UAA0 + ZM1*UAA1 + ZM2*UAA2 + ZM3*UAA3) / (tet->Det()*6.0)
00627           * Coefficient_triangle;
00628       }
00629       else if(tet->N1()==num_v) {
00630         sum = sum +
00631           ZM1*(ZM0*UAA0 + ZM1*UAA1 + ZM2*UAA2 + ZM3*UAA3) / (tet->Det()*6.0)
00632           * Coefficient_triangle;
00633       }
00634       else if(tet->N2()==num_v) {
00635         sum = sum +
00636           ZM2*(ZM0*UAA0 + ZM1*UAA1 + ZM2*UAA2 + ZM3*UAA3) / (tet->Det()*6.0)
00637           * Coefficient_triangle;
00638       }
00639       else if(tet->N3()==num_v) {
00640         sum = sum +
00641           ZM3*(ZM0*UAA0 + ZM1*UAA1 + ZM2*UAA2 + ZM3*UAA3) / (tet->Det()*6.0)
00642           * Coefficient_triangle;
00643       }
00644     }     
00645     return sum;
00646   }
00647   static inline double diag_F_bo_cell(const BoCeData* bocedata,
00648                                       int num_v, int num_coefficient) {
00649    static double sum;
00650     Tetraeder_storage* tet;
00651 
00652     sum=0.0;
00653     for(tet=bocedata->Give_tets();tet!=NULL;tet=tet->Next()) {
00654       if(tet->N0()==num_v) {
00655         sum = sum + ZM0*ZM0 / (tet->Det()*6.0)
00656           * Coefficient_triangle;
00657       }
00658       else if(tet->N1()==num_v) {
00659         sum = sum + ZM1*ZM1 / (tet->Det()*6.0)
00660           * Coefficient_triangle;
00661       }
00662       else if(tet->N2()==num_v) {
00663         sum = sum + ZM2*ZM2 / (tet->Det()*6.0)
00664           * Coefficient_triangle;
00665       }
00666       else if(tet->N3()==num_v) {
00667         sum = sum + ZM3*ZM3 / (tet->Det()*6.0)
00668           * Coefficient_triangle;
00669       }
00670     }     
00671     return sum;
00672   }
00673 };
00674 
00675 
00676 //  Helm
00677 class helm_FE_variable {
00678  public:
00679   diff_type typ_u() { return  Helm_type; }; 
00680   diff_type typ_v() { return  Helm_type; };
00681 
00682   helm_FE_variable() { }
00683   static inline int    ops_interior() { return 88; }
00684   static inline int    ops_diag_interior() { return 21; }
00685   static inline double stencil(double uN,
00686                                double uW, double uM, double uE,
00687                                double uS, double uT, double uD,
00688                                double uND, double uWN, double uWT,
00689                                double uED, double uST, double uES,
00690                                double uEST, double uWND,
00691                                double aN,
00692                                double aW, double aM, double aE,
00693                                double aS, double aT, double aD,
00694                                double aST, double aET, double aES,
00695                                double aND, double aWD, double aWN,
00696                                double aEST, double aWND,
00697                                double h_mesh) {
00698    return (
00699            var_interpo_WSD * (uW + 2.0*uM + uD + uS) +
00700            var_interpo_ESD * (6.0*uM + 2.0*(uED + uD) + uE + uS + 3.0*uES) +
00701            var_interpo_WND * (3.0*uWN + 4.0*uWND + uN + 3.0*uND + 
00702                               10.0*uM + 2.0*(uD + uW)) +
00703            var_interpo_END * (uN + 3.0*uND + 6.0*uM + uD + 2.0*(uE + uED)) +
00704            var_interpo_WST * (2.0*(uWT + uW) + uT + 6.0*uM + 3.0*uST + uS) +
00705            var_interpo_EST * (2.0*(uT + uE) + 3.0*(uST + uES) + uS + 
00706                               4.0*uEST + 10.0*uM) +
00707            var_interpo_WNT * (3.0*uWN + uN + 2.0*(uWT + uT) + uW + 6.0*uM) +
00708            var_interpo_ENT * (uN + uT + 2.0*uM + uE)
00709            ) * h_mesh*h_mesh*h_mesh / 240.0; // additional divide by 2.0
00710   }
00711   static inline double diag_stencil(double h_mesh,
00712                                     double aN,
00713                                     double aW, double aM, double aE,
00714                                     double aS, double aT, double aD,
00715                                     double aST, double aET, double aES,
00716                                     double aND, double aWD, double aWN,
00717                                     double aEST, double aWND) {
00718    return (
00719            (var_interpo_EST +
00720             var_interpo_WND) * 5.0 + 
00721            (var_interpo_ESD +
00722             var_interpo_END +
00723             var_interpo_WST +
00724             var_interpo_WNT) * 3.0 +
00725            (var_interpo_ENT  +
00726             var_interpo_WSD)
00727            ) * h_mesh*h_mesh*h_mesh / 120.0; // additional divide by 2.0
00728   }
00729   static inline double Factor_reg_cell(double h_mesh) {
00730     return h_mesh*h_mesh*h_mesh / 120.0;
00731   }
00732   static inline double diag_F_reg_cell(dir_sons dir_v,
00733                                        double h_mesh, double* sten) {
00734     return sten[dir_v*9] * h_mesh*h_mesh*h_mesh / 120.0;;
00735   }
00736 
00737   // Operator applied to variable with number
00738   static inline double F_bo_cell(const BoCeData* bocedata,
00739                                  int num_v, int num_variable,
00740                                  int num_coefficient) {
00741     static double sum;
00742     Tetraeder_storage* tet;
00743     
00744     sum=0.0;
00745     for(tet=bocedata->Give_tets();tet!=NULL;tet=tet->Next()) {
00746       if(tet->N0()==num_v) {
00747         sum = sum + tet->Det() *
00748               (2.0 * UU0 + UU1 + UU2 + UU3) / 120.0
00749           * Coefficient_triangle;
00750       }
00751       else if(tet->N1()==num_v) {
00752         sum = sum + tet->Det() *
00753               (2.0 * UU1 + UU0 + UU2 + UU3) / 120.0
00754           * Coefficient_triangle;
00755       }
00756       else if(tet->N2()==num_v) {
00757         sum = sum + tet->Det() *
00758               (2.0 * UU2 + UU0 + UU1 + UU3) / 120.0
00759           * Coefficient_triangle;
00760       }
00761       else if(tet->N3()==num_v) {
00762         sum = sum + tet->Det() *
00763               (2.0 * UU3 + UU0 + UU1 + UU2) / 120.0
00764           * Coefficient_triangle;
00765       }
00766     }  
00767     return sum;
00768   }
00769   // Operator applied to a vector U_array
00770   static inline double Local_matrix_bo_cell(const BoCeData* bocedata,
00771                                             double*  U_array, int num_v,
00772                                             int num_coefficient) {
00773     static double sum;
00774     Tetraeder_storage* tet;
00775     
00776     sum=0.0;
00777     for(tet=bocedata->Give_tets();tet!=NULL;tet=tet->Next()) {
00778       if(tet->N0()==num_v) {
00779         sum = sum + tet->Det() *
00780               (2.0 * UAA0 + UAA1 + UAA2 + UAA3) / 120.0
00781           * Coefficient_triangle;
00782       }
00783       else if(tet->N1()==num_v) {
00784         sum = sum + tet->Det() *
00785               (2.0 * UAA1 + UAA0 + UAA2 + UAA3) / 120.0
00786           * Coefficient_triangle;
00787       }
00788       else if(tet->N2()==num_v) {
00789         sum = sum + tet->Det() *
00790               (2.0 * UAA2 + UAA0 + UAA1 + UAA3) / 120.0
00791           * Coefficient_triangle;
00792       }
00793       else if(tet->N3()==num_v) {
00794         sum = sum + tet->Det() *
00795               (2.0 * UAA3 + UAA0 + UAA1 + UAA2) / 120.0
00796           * Coefficient_triangle;
00797       }
00798     }  
00799     return sum;
00800   }
00801   static inline double diag_F_bo_cell(const BoCeData* bocedata,
00802                                       int num_v, int num_coefficient) {
00803    static double sum;
00804     Tetraeder_storage* tet;
00805     
00806     sum=0.0;
00807     for(tet=bocedata->Give_tets();tet!=NULL;tet=tet->Next()) {
00808       if(tet->N0()==num_v || tet->N1()==num_v ||
00809          tet->N2()==num_v || tet->N3()==num_v)
00810         sum = sum + tet->Det() / 60.0
00811           * Coefficient_triangle;
00812     }  
00813     return sum;
00814   }
00815 };
00816 
00818 // Unsymmetrische Terme
00820 
00821 
00822 // Formeln fuer dxdy
00823 class dxdy_FE_variable {
00824  public:  
00825   diff_type typ_u() { return  Dx_type; }; 
00826   diff_type typ_v() { return  Dy_type; };
00827 
00828   dxdy_FE_variable() { }
00829   static inline int    ops_interior() { return 46; } // 38 + 8
00830   static inline int    ops_diag_interior() { return 16; }
00831   static inline double stencil(double uN,
00832                                double uW, double uM, double uE,
00833                                double uS, double uT, double uD,
00834                                double uND, double uWN, double uWT,
00835                                double uED, double uST, double uES,
00836                                double uEST, double uWND,
00837                                double aN,
00838                                double aW, double aM, double aE,
00839                                double aS, double aT, double aD,
00840                                double aST, double aET, double aES,
00841                                double aND, double aWD, double aWN,
00842                                double aEST, double aWND,
00843                                double h_mesh) {
00844     return (
00845             var_interpo_WSD * ( uM - uW) +
00846             var_interpo_ESD * ( uED - uD - uS + uES) +
00847             var_interpo_WND * (   uWN + uWND - uN - uND - uW + uM) +
00848             (var_interpo_END * (   uM - uE) +
00849              var_interpo_WST * ( uM - uW) ) * 2.0 +
00850             var_interpo_EST * (   uM - uE - uST - uS + uEST + uES) +
00851             var_interpo_WNT * (   uWN - uN + uWT - uT) +
00852             var_interpo_ENT * (   uM - uE)
00853             ) * h_mesh / 12.0; 
00854   }
00855   static inline double diag_stencil(double h_mesh,
00856                                     double aN,
00857                                     double aW, double aM, double aE,
00858                                     double aS, double aT, double aD,
00859                                     double aST, double aET, double aES,
00860                                     double aND, double aWD, double aWN,
00861                                     double aEST, double aWND) {
00862     return (
00863             var_interpo_WSD +
00864             var_interpo_WND +
00865             (var_interpo_END +
00866              var_interpo_WST) * 2.0 +
00867             var_interpo_EST +
00868             var_interpo_ENT 
00869             ) * h_mesh / 12.0; 
00870   }
00871   static inline double Factor_reg_cell(double h_mesh) {
00872     return h_mesh / 6.0;
00873   }
00874   static inline double diag_F_reg_cell(dir_sons dir_v,
00875                                        double h_mesh, double* sten) {
00876     return sten[dir_v*9] * h_mesh / 6.0;
00877   }
00878 
00879   // Operator applied to variable with number
00880   static inline double F_bo_cell(const BoCeData* bocedata,
00881                                  int num_v, int num_variable,
00882                                  int num_coefficient) {
00883     static double sum;
00884     Tetraeder_storage* tet;
00885 
00886     sum=0.0;
00887     for(tet=bocedata->Give_tets();tet!=NULL;tet=tet->Next()) {
00888       if(tet->N0()==num_v) {
00889         sum = sum +
00890               YM0*(XM0*UU0 + XM1*UU1 + XM2*UU2 + XM3*UU3) / (tet->Det()*6.0)
00891           * Coefficient_triangle;
00892       }
00893       else if(tet->N1()==num_v) {
00894         sum = sum +
00895               YM1*(XM0*UU0 + XM1*UU1 + XM2*UU2 + XM3*UU3) / (tet->Det()*6.0)
00896           * Coefficient_triangle;
00897       }
00898       else if(tet->N2()==num_v) {
00899         sum = sum +
00900               YM2*(XM0*UU0 + XM1*UU1 + XM2*UU2 + XM3*UU3) / (tet->Det()*6.0)
00901           * Coefficient_triangle;
00902       }
00903       else if(tet->N3()==num_v) {
00904         sum = sum +
00905               YM3*(XM0*UU0 + XM1*UU1 + XM2*UU2 + XM3*UU3) / (tet->Det()*6.0)
00906           * Coefficient_triangle;
00907       }
00908     }     
00909     return sum;
00910   }
00911   // Operator applied to a vector U_array
00912   static inline double Local_matrix_bo_cell(const BoCeData* bocedata,
00913                                             double*  U_array, int num_v,
00914                                             int num_coefficient) {
00915     static double sum;
00916     Tetraeder_storage* tet;
00917 
00918     sum=0.0;
00919     for(tet=bocedata->Give_tets();tet!=NULL;tet=tet->Next()) {
00920       if(tet->N0()==num_v) {
00921         sum = sum +
00922           YM0*(XM0*UAA0 + XM1*UAA1 + XM2*UAA2 + XM3*UAA3) / (tet->Det()*6.0)
00923           * Coefficient_triangle;
00924       }
00925       else if(tet->N1()==num_v) {
00926         sum = sum +
00927           YM1*(XM0*UAA0 + XM1*UAA1 + XM2*UAA2 + XM3*UAA3) / (tet->Det()*6.0)
00928           * Coefficient_triangle;
00929       }
00930       else if(tet->N2()==num_v) {
00931         sum = sum +
00932           YM2*(XM0*UAA0 + XM1*UAA1 + XM2*UAA2 + XM3*UAA3) / (tet->Det()*6.0)
00933           * Coefficient_triangle;
00934       }
00935       else if(tet->N3()==num_v) {
00936         sum = sum +
00937           YM3*(XM0*UAA0 + XM1*UAA1 + XM2*UAA2 + XM3*UAA3) / (tet->Det()*6.0)
00938           * Coefficient_triangle;
00939       }
00940     }     
00941     return sum;
00942   }
00943   static inline double diag_F_bo_cell(const BoCeData* bocedata,
00944                                       int num_v, int num_coefficient) {
00945    static double sum;
00946     Tetraeder_storage* tet;
00947 
00948     sum=0.0;
00949     for(tet=bocedata->Give_tets();tet!=NULL;tet=tet->Next()) {
00950       if(tet->N0()==num_v) {
00951         sum = sum + YM0*XM0 / (tet->Det()*6.0)
00952           * Coefficient_triangle;
00953       }
00954       else if(tet->N1()==num_v) {
00955         sum = sum + YM1*XM1 / (tet->Det()*6.0)
00956           * Coefficient_triangle;
00957       }
00958       else if(tet->N2()==num_v) {
00959         sum = sum + YM2*XM2 / (tet->Det()*6.0)
00960           * Coefficient_triangle;
00961       }
00962       else if(tet->N3()==num_v) {
00963         sum = sum + YM3*XM3 / (tet->Det()*6.0)
00964           * Coefficient_triangle;
00965       }
00966     }     
00967     return sum;
00968   }
00969 };
00970 
00971 // Formeln fuer dydx
00972 class dydx_FE_variable {
00973  public:  
00974   diff_type typ_u() { return  Dy_type; }; 
00975   diff_type typ_v() { return  Dx_type; };
00976 
00977   dydx_FE_variable() { }
00978   static inline int    ops_interior() { return 41; }
00979   static inline int    ops_diag_interior() { return 16; }
00980   static inline double stencil(double uN,
00981                                double uW, double uM, double uE,
00982                                double uS, double uT, double uD,
00983                                double uND, double uWN, double uWT,
00984                                double uED, double uST, double uES,
00985                                double uEST, double uWND,
00986                                double aN,
00987                                double aW, double aM, double aE,
00988                                double aS, double aT, double aD,
00989                                double aST, double aET, double aES,
00990                                double aND, double aWD, double aWN,
00991                                double aEST, double aWND,
00992                                double h_mesh) {
00993     return (
00994             var_interpo_WSD * (   uM - uS) +
00995             var_interpo_ESD * (   uES - uE) +
00996             var_interpo_WND * (   uWN + uWND - 2.0*uW + uM - uD) +
00997             var_interpo_END * ( 2.0*uM - uN - uND - uE + uED) +
00998             var_interpo_WST * (   uWT - uW + 2.0*uM - uST - uS) +
00999             var_interpo_EST * ( uM - uT - 2.0*uE + uEST + uES) +
01000             var_interpo_WNT * (   uWN - uW) +
01001             var_interpo_ENT * ( uM - uN)
01002             )  * h_mesh / 12.0; 
01003   }
01004   static inline double diag_stencil(double h_mesh,
01005                                     double aN,
01006                                     double aW, double aM, double aE,
01007                                     double aS, double aT, double aD,
01008                                     double aST, double aET, double aES,
01009                                     double aND, double aWD, double aWN,
01010                                     double aEST, double aWND) {
01011     return (
01012             var_interpo_WSD +
01013             var_interpo_WND +
01014             (var_interpo_END +
01015              var_interpo_WST) * 2.0 +
01016             var_interpo_EST +
01017             var_interpo_ENT
01018             )  * h_mesh / 12.0; 
01019   }
01020   static inline double Factor_reg_cell(double h_mesh) {
01021     return h_mesh / 6.0;
01022   }
01023   static inline double diag_F_reg_cell(dir_sons dir_v,
01024                                        double h_mesh, double* sten) {
01025     return sten[dir_v*9] * h_mesh / 6.0;
01026   }
01027 
01028   // Operator applied to variable with number
01029   static inline double F_bo_cell(const BoCeData* bocedata,
01030                                  int num_v, int num_variable,
01031                                  int num_coefficient) {
01032     static double sum;
01033     Tetraeder_storage* tet;
01034 
01035     sum=0.0;
01036     for(tet=bocedata->Give_tets();tet!=NULL;tet=tet->Next()) {
01037       if(tet->N0()==num_v) {
01038         sum = sum +
01039               XM0*(YM0*UU0 + YM1*UU1 + YM2*UU2 + YM3*UU3) / (tet->Det()*6.0)
01040           * Coefficient_triangle;
01041       }
01042       else if(tet->N1()==num_v) {
01043         sum = sum +
01044               XM1*(YM0*UU0 + YM1*UU1 + YM2*UU2 + YM3*UU3) / (tet->Det()*6.0)
01045           * Coefficient_triangle;
01046       }
01047       else if(tet->N2()==num_v) {
01048         sum = sum +
01049               XM2*(YM0*UU0 + YM1*UU1 + YM2*UU2 + YM3*UU3) / (tet->Det()*6.0)
01050           * Coefficient_triangle;
01051       }
01052       else if(tet->N3()==num_v) {
01053         sum = sum +
01054               XM3*(YM0*UU0 + YM1*UU1 + YM2*UU2 + YM3*UU3) / (tet->Det()*6.0)
01055           * Coefficient_triangle;
01056       }
01057     }     
01058     return sum;
01059   }
01060   // Operator applied to a vector U_array
01061   static inline double Local_matrix_bo_cell(const BoCeData* bocedata,
01062                                             double*  U_array, int num_v,
01063                                             int num_coefficient) {
01064     static double sum;
01065     Tetraeder_storage* tet;
01066 
01067     sum=0.0;
01068     for(tet=bocedata->Give_tets();tet!=NULL;tet=tet->Next()) {
01069       if(tet->N0()==num_v) {
01070         sum = sum +
01071           XM0*(YM0*UAA0 + YM1*UAA1 + YM2*UAA2 + YM3*UAA3) / (tet->Det()*6.0)
01072           * Coefficient_triangle;
01073       }
01074       else if(tet->N1()==num_v) {
01075         sum = sum +
01076           XM1*(YM0*UAA0 + YM1*UAA1 + YM2*UAA2 + YM3*UAA3) / (tet->Det()*6.0)
01077           * Coefficient_triangle;
01078       }
01079       else if(tet->N2()==num_v) {
01080         sum = sum +
01081           XM2*(YM0*UAA0 + YM1*UAA1 + YM2*UAA2 + YM3*UAA3) / (tet->Det()*6.0)
01082           * Coefficient_triangle;
01083       }
01084       else if(tet->N3()==num_v) {
01085         sum = sum +
01086           XM3*(YM0*UAA0 + YM1*UAA1 + YM2*UAA2 + YM3*UAA3) / (tet->Det()*6.0)
01087           * Coefficient_triangle;
01088       }
01089     }     
01090     return sum;
01091   }
01092   static inline double diag_F_bo_cell(const BoCeData* bocedata,
01093                                       int num_v, int num_coefficient) {
01094    static double sum;
01095     Tetraeder_storage* tet;
01096 
01097     sum=0.0;
01098     for(tet=bocedata->Give_tets();tet!=NULL;tet=tet->Next()) {
01099       if(tet->N0()==num_v) {
01100         sum = sum + XM0*YM0 / (tet->Det()*6.0)
01101           * Coefficient_triangle;
01102       }
01103       else if(tet->N1()==num_v) {
01104         sum = sum + XM1*YM1 / (tet->Det()*6.0)
01105           * Coefficient_triangle;
01106       }
01107       else if(tet->N2()==num_v) {
01108         sum = sum + XM2*YM2 / (tet->Det()*6.0)
01109           * Coefficient_triangle;
01110       }
01111       else if(tet->N3()==num_v) {
01112         sum = sum + XM3*YM3 / (tet->Det()*6.0)
01113           * Coefficient_triangle;
01114       }
01115     }     
01116     return sum;
01117   }
01118 };
01119 
01121 
01122 
01123 
01124 // Formeln fuer dxdz
01125 class dxdz_FE_variable {
01126  public:  
01127   diff_type typ_u() { return  Dx_type; }; 
01128   diff_type typ_v() { return  Dz_type; };
01129 
01130   dxdz_FE_variable() { }
01131   static inline int    ops_interior() { return 41; } // 33 + 8
01132   static inline int    ops_diag_interior() { return 13; }
01133   static inline double stencil(double uN,
01134                                double uW, double uM, double uE,
01135                                double uS, double uT, double uD,
01136                                double uND, double uWN, double uWT,
01137                                double uED, double uST, double uES,
01138                                double uEST, double uWND,
01139                                double aN,
01140                                double aW, double aM, double aE,
01141                                double aS, double aT, double aD,
01142                                double aST, double aET, double aES,
01143                                double aND, double aWD, double aWN,
01144                                double aEST, double aWND,
01145                                double h_mesh) {
01146     return (
01147             var_interpo_WSD * ( uM - uW) +
01148             var_interpo_ESD * ( uED - uD - uS + uES) +
01149             var_interpo_WND * ( uM - uWND + uND - uW) +
01150             var_interpo_END * ( uED - uD) +
01151             var_interpo_WST * (   uWT - uT) +
01152             var_interpo_EST * (   uM - uE + uST - uEST) +
01153             var_interpo_WNT * (   uWN - uN + uWT - uT) +
01154             var_interpo_ENT * (   uM - uE)
01155             ) * h_mesh / 12.0; 
01156   }
01157   static inline double diag_stencil(double h_mesh,
01158                                     double aN,
01159                                     double aW, double aM, double aE,
01160                                     double aS, double aT, double aD,
01161                                     double aST, double aET, double aES,
01162                                     double aND, double aWD, double aWN,
01163                                     double aEST, double aWND) {
01164     return (
01165             var_interpo_WSD +
01166             var_interpo_WND +
01167             var_interpo_EST +
01168             var_interpo_ENT
01169             ) * h_mesh / 12.0; 
01170   }
01171   static inline double Factor_reg_cell(double h_mesh) {
01172     return h_mesh / 6.0;
01173   }
01174   static inline double diag_F_reg_cell(dir_sons dir_v,
01175                                        double h_mesh, double* sten) {
01176     return sten[dir_v*9] * h_mesh / 6.0;
01177   }
01178 
01179   // Operator applied to variable with number
01180   static inline double F_bo_cell(const BoCeData* bocedata,
01181                                  int num_v, int num_variable,
01182                                  int num_coefficient) {
01183     static double sum;
01184     Tetraeder_storage* tet;
01185 
01186     sum=0.0;
01187     for(tet=bocedata->Give_tets();tet!=NULL;tet=tet->Next()) {
01188       if(tet->N0()==num_v) {
01189         sum = sum +
01190               ZM0*(XM0*UU0 + XM1*UU1 + XM2*UU2 + XM3*UU3) / (tet->Det()*6.0)
01191           * Coefficient_triangle;
01192       }
01193       else if(tet->N1()==num_v) {
01194         sum = sum +
01195               ZM1*(XM0*UU0 + XM1*UU1 + XM2*UU2 + XM3*UU3) / (tet->Det()*6.0)
01196           * Coefficient_triangle;
01197       }
01198       else if(tet->N2()==num_v) {
01199         sum = sum +
01200               ZM2*(XM0*UU0 + XM1*UU1 + XM2*UU2 + XM3*UU3) / (tet->Det()*6.0)
01201           * Coefficient_triangle;
01202       }
01203       else if(tet->N3()==num_v) {
01204         sum = sum +
01205               ZM3*(XM0*UU0 + XM1*UU1 + XM2*UU2 + XM3*UU3) / (tet->Det()*6.0)
01206           * Coefficient_triangle;
01207       }
01208     }     
01209     return sum;
01210   }
01211   // Operator applied to a vector U_array
01212   static inline double Local_matrix_bo_cell(const BoCeData* bocedata,
01213                                             double*  U_array, int num_v,
01214                                             int num_coefficient) {
01215     static double sum;
01216     Tetraeder_storage* tet;
01217 
01218     sum=0.0;
01219     for(tet=bocedata->Give_tets();tet!=NULL;tet=tet->Next()) {
01220       if(tet->N0()==num_v) {
01221         sum = sum +
01222           ZM0*(XM0*UAA0 + XM1*UAA1 + XM2*UAA2 + XM3*UAA3) / (tet->Det()*6.0)
01223           * Coefficient_triangle;
01224       }
01225       else if(tet->N1()==num_v) {
01226         sum = sum +
01227           ZM1*(XM0*UAA0 + XM1*UAA1 + XM2*UAA2 + XM3*UAA3) / (tet->Det()*6.0)
01228           * Coefficient_triangle;
01229       }
01230       else if(tet->N2()==num_v) {
01231         sum = sum +
01232           ZM2*(XM0*UAA0 + XM1*UAA1 + XM2*UAA2 + XM3*UAA3) / (tet->Det()*6.0)
01233           * Coefficient_triangle;
01234       }
01235       else if(tet->N3()==num_v) {
01236         sum = sum +
01237           ZM3*(XM0*UAA0 + XM1*UAA1 + XM2*UAA2 + XM3*UAA3) / (tet->Det()*6.0)
01238           * Coefficient_triangle;
01239       }
01240     }     
01241     return sum;
01242   }
01243   static inline double diag_F_bo_cell(const BoCeData* bocedata,
01244                                       int num_v, int num_coefficient) {
01245    static double sum;
01246     Tetraeder_storage* tet;
01247 
01248     sum=0.0;
01249     for(tet=bocedata->Give_tets();tet!=NULL;tet=tet->Next()) {
01250       if(tet->N0()==num_v) {
01251         sum = sum + ZM0*XM0 / (tet->Det()*6.0)
01252           * Coefficient_triangle;
01253       }
01254       else if(tet->N1()==num_v) {
01255         sum = sum + ZM1*XM1 / (tet->Det()*6.0)
01256           * Coefficient_triangle;
01257       }
01258       else if(tet->N2()==num_v) {
01259         sum = sum + ZM2*XM2 / (tet->Det()*6.0)
01260           * Coefficient_triangle;
01261       }
01262       else if(tet->N3()==num_v) {
01263         sum = sum + ZM3*XM3 / (tet->Det()*6.0)
01264           * Coefficient_triangle;
01265       }
01266     }     
01267     return sum;
01268   }
01269 };
01270 
01271 // Formeln fuer dzdx
01272 class dzdx_FE_variable {
01273  public:  
01274   diff_type typ_u() { return  Dz_type; }; 
01275   diff_type typ_v() { return  Dx_type; };
01276 
01277   dzdx_FE_variable() { }
01278   static inline int    ops_interior() { return 41; }  // 33 + 8
01279   static inline int    ops_diag_interior() { return 13; }
01280   static inline double stencil(double uN,
01281                                double uW, double uM, double uE,
01282                                double uS, double uT, double uD,
01283                                double uND, double uWN, double uWT,
01284                                double uED, double uST, double uES,
01285                                double uEST, double uWND,
01286                                double aN,
01287                                double aW, double aM, double aE,
01288                                double aS, double aT, double aD,
01289                                double aST, double aET, double aES,
01290                                double aND, double aWD, double aWN,
01291                                double aEST, double aWND,
01292                                double h_mesh) {
01293     return (
01294             var_interpo_WSD * (   uM - uD) +
01295             var_interpo_ESD * ( - uE + uED) +
01296             var_interpo_WND * (   uWN - uWND + uM - uD) +
01297             var_interpo_END * ( - uN + uND - uE + uED) +
01298             var_interpo_WST * (   uWT - uW + uST - uS) +
01299             var_interpo_EST * ( - uT + uM - uEST + uES) +
01300             var_interpo_WNT * (   uWT - uW) +
01301             var_interpo_ENT * ( - uT + uM)
01302             ) * h_mesh / 12.0;   
01303   }
01304   static inline double diag_stencil(double h_mesh,
01305                                     double aN,
01306                                     double aW, double aM, double aE,
01307                                     double aS, double aT, double aD,
01308                                     double aST, double aET, double aES,
01309                                     double aND, double aWD, double aWN,
01310                                     double aEST, double aWND) {
01311     return (
01312             var_interpo_WSD +
01313             var_interpo_WND +
01314             var_interpo_EST +
01315             var_interpo_ENT 
01316             ) * h_mesh / 12.0;   
01317   }
01318   static inline double Factor_reg_cell(double h_mesh) {
01319     return h_mesh / 6.0;
01320   }
01321   static inline double diag_F_reg_cell(dir_sons dir_v,
01322                                        double h_mesh, double* sten) {
01323     return sten[dir_v*9] * h_mesh / 6.0;
01324   }
01325 
01326   // Operator applied to variable with number
01327   static inline double F_bo_cell(const BoCeData* bocedata,
01328                                  int num_v, int num_variable,
01329                                  int num_coefficient) {
01330     static double sum;
01331     Tetraeder_storage* tet;
01332 
01333     sum=0.0;
01334     for(tet=bocedata->Give_tets();tet!=NULL;tet=tet->Next()) {
01335       if(tet->N0()==num_v) {
01336         sum = sum +
01337               XM0*(ZM0*UU0 + ZM1*UU1 + ZM2*UU2 + ZM3*UU3) / (tet->Det()*6.0)
01338           * Coefficient_triangle;
01339       }
01340       else if(tet->N1()==num_v) {
01341         sum = sum +
01342               XM1*(ZM0*UU0 + ZM1*UU1 + ZM2*UU2 + ZM3*UU3) / (tet->Det()*6.0)
01343           * Coefficient_triangle;
01344       }
01345       else if(tet->N2()==num_v) {
01346         sum = sum +
01347               XM2*(ZM0*UU0 + ZM1*UU1 + ZM2*UU2 + ZM3*UU3) / (tet->Det()*6.0)
01348           * Coefficient_triangle;
01349       }
01350       else if(tet->N3()==num_v) {
01351         sum = sum +
01352               XM3*(ZM0*UU0 + ZM1*UU1 + ZM2*UU2 + ZM3*UU3) / (tet->Det()*6.0)
01353           * Coefficient_triangle;
01354       }
01355     }     
01356     return sum;
01357   }
01358   // Operator applied to a vector U_array
01359   static inline double Local_matrix_bo_cell(const BoCeData* bocedata,
01360                                             double*  U_array, int num_v,
01361                                             int num_coefficient) {
01362     static double sum;
01363     Tetraeder_storage* tet;
01364 
01365     sum=0.0;
01366     for(tet=bocedata->Give_tets();tet!=NULL;tet=tet->Next()) {
01367       if(tet->N0()==num_v) {
01368         sum = sum +
01369           XM0*(ZM0*UAA0 + ZM1*UAA1 + ZM2*UAA2 + ZM3*UAA3) / (tet->Det()*6.0)
01370           * Coefficient_triangle;
01371       }
01372       else if(tet->N1()==num_v) {
01373         sum = sum +
01374           XM1*(ZM0*UAA0 + ZM1*UAA1 + ZM2*UAA2 + ZM3*UAA3) / (tet->Det()*6.0)
01375           * Coefficient_triangle;
01376       }
01377       else if(tet->N2()==num_v) {
01378         sum = sum +
01379           XM2*(ZM0*UAA0 + ZM1*UAA1 + ZM2*UAA2 + ZM3*UAA3) / (tet->Det()*6.0)
01380           * Coefficient_triangle;
01381       }
01382       else if(tet->N3()==num_v) {
01383         sum = sum +
01384           XM3*(ZM0*UAA0 + ZM1*UAA1 + ZM2*UAA2 + ZM3*UAA3) / (tet->Det()*6.0)
01385           * Coefficient_triangle;
01386       }
01387     }     
01388     return sum;
01389   }
01390   static inline double diag_F_bo_cell(const BoCeData* bocedata,
01391                                       int num_v, int num_coefficient) {
01392    static double sum;
01393     Tetraeder_storage* tet;
01394 
01395     sum=0.0;
01396     for(tet=bocedata->Give_tets();tet!=NULL;tet=tet->Next()) {
01397       if(tet->N0()==num_v) {
01398         sum = sum + XM0*ZM0 / (tet->Det()*6.0)
01399           * Coefficient_triangle;
01400       }
01401       else if(tet->N1()==num_v) {
01402         sum = sum + XM1*ZM1 / (tet->Det()*6.0)
01403           * Coefficient_triangle;
01404       }
01405       else if(tet->N2()==num_v) {
01406         sum = sum + XM2*ZM2 / (tet->Det()*6.0)
01407           * Coefficient_triangle;
01408       }
01409       else if(tet->N3()==num_v) {
01410         sum = sum + XM3*ZM3 / (tet->Det()*6.0)
01411           * Coefficient_triangle;
01412       }
01413     }     
01414     return sum;
01415   }
01416 };
01417 
01418 
01420 
01421 
01422 // Formeln fuer dydz
01423 class dydz_FE_variable {
01424  public:  
01425   diff_type typ_u() { return  Dy_type; }; 
01426   diff_type typ_v() { return  Dz_type; };
01427 
01428   dydz_FE_variable() { }
01429   static inline int    ops_interior() { return 49; } // 41 + 8
01430   static inline int    ops_diag_interior() { return 16; }
01431   static inline double stencil(double uN,
01432                                double uW, double uM, double uE,
01433                                double uS, double uT, double uD,
01434                                double uND, double uWN, double uWT,
01435                                double uED, double uST, double uES,
01436                                double uEST, double uWND,
01437                                double aN,
01438                                double aW, double aM, double aE,
01439                                double aS, double aT, double aD,
01440                                double aST, double aET, double aES,
01441                                double aND, double aWD, double aWN,
01442                                double aEST, double aWND,
01443                                double h_mesh) {
01444     return (
01445             var_interpo_WSD * (   uM - uS) +
01446             var_interpo_ESD * (   2.0*uM - uD + uED - uS - uES) +
01447             var_interpo_WND * (   uWND + uND - uW + uM - 2.0*uD) +
01448             var_interpo_END * (   uND - uD) +
01449             var_interpo_WST * ( - uT + uST) +
01450             var_interpo_EST * ( - 2.0*uT + uM - uE + uST + uEST) +
01451             var_interpo_WNT * ( - uWN - uN + uWT - uT + 2.0*uM) +
01452             var_interpo_ENT * ( - uN + uM)
01453             ) * h_mesh / 12.0; 
01454   }
01455   static inline double diag_stencil(double h_mesh,
01456                                     double aN,
01457                                     double aW, double aM, double aE,
01458                                     double aS, double aT, double aD,
01459                                     double aST, double aET, double aES,
01460                                     double aND, double aWD, double aWN,
01461                                     double aEST, double aWND) {
01462     return (
01463             var_interpo_WSD +
01464             (var_interpo_ESD +
01465              var_interpo_WNT) * 2.0 +
01466             var_interpo_WND +
01467             var_interpo_EST +
01468             var_interpo_ENT 
01469             ) * h_mesh / 12.0; 
01470   }
01471   static inline double Factor_reg_cell(double h_mesh) {
01472     return h_mesh / 6.0;
01473   }
01474   static inline double diag_F_reg_cell(dir_sons dir_v,
01475                                        double h_mesh, double* sten) {
01476     return sten[dir_v*9] * h_mesh / 6.0;
01477   }
01478 
01479   // Operator applied to variable with number
01480   static inline double F_bo_cell(const BoCeData* bocedata,
01481                                  int num_v, int num_variable,
01482                                  int num_coefficient) {
01483     static double sum;
01484     Tetraeder_storage* tet;
01485 
01486     sum=0.0;
01487     for(tet=bocedata->Give_tets();tet!=NULL;tet=tet->Next()) {
01488       if(tet->N0()==num_v) {
01489         sum = sum +
01490               ZM0*(YM0*UU0 + YM1*UU1 + YM2*UU2 + YM3*UU3) / (tet->Det()*6.0)
01491           * Coefficient_triangle;
01492       }
01493       else if(tet->N1()==num_v) {
01494         sum = sum +
01495               ZM1*(YM0*UU0 + YM1*UU1 + YM2*UU2 + YM3*UU3) / (tet->Det()*6.0)
01496           * Coefficient_triangle;
01497       }
01498       else if(tet->N2()==num_v) {
01499         sum = sum +
01500               ZM2*(YM0*UU0 + YM1*UU1 + YM2*UU2 + YM3*UU3) / (tet->Det()*6.0)
01501           * Coefficient_triangle;
01502       }
01503       else if(tet->N3()==num_v) {
01504         sum = sum +
01505               ZM3*(YM0*UU0 + YM1*UU1 + YM2*UU2 + YM3*UU3) / (tet->Det()*6.0)
01506           * Coefficient_triangle;
01507       }
01508     }     
01509     return sum;
01510   }
01511   // Operator applied to a vector U_array
01512   static inline double Local_matrix_bo_cell(const BoCeData* bocedata,
01513                                             double*  U_array, int num_v,
01514                                             int num_coefficient) {
01515     static double sum;
01516     Tetraeder_storage* tet;
01517 
01518     sum=0.0;
01519     for(tet=bocedata->Give_tets();tet!=NULL;tet=tet->Next()) {
01520       if(tet->N0()==num_v) {
01521         sum = sum +
01522           ZM0*(YM0*UAA0 + YM1*UAA1 + YM2*UAA2 + YM3*UAA3) / (tet->Det()*6.0)
01523           * Coefficient_triangle;
01524       }
01525       else if(tet->N1()==num_v) {
01526         sum = sum +
01527           ZM1*(YM0*UAA0 + YM1*UAA1 + YM2*UAA2 + YM3*UAA3) / (tet->Det()*6.0)
01528           * Coefficient_triangle;
01529       }
01530       else if(tet->N2()==num_v) {
01531         sum = sum +
01532           ZM2*(YM0*UAA0 + YM1*UAA1 + YM2*UAA2 + YM3*UAA3) / (tet->Det()*6.0)
01533           * Coefficient_triangle;
01534       }
01535       else if(tet->N3()==num_v) {
01536         sum = sum +
01537           ZM3*(YM0*UAA0 + YM1*UAA1 + YM2*UAA2 + YM3*UAA3) / (tet->Det()*6.0)
01538           * Coefficient_triangle;
01539       }
01540     }     
01541     return sum;
01542   }
01543   static inline double diag_F_bo_cell(const BoCeData* bocedata,
01544                                       int num_v, int num_coefficient) {
01545    static double sum;
01546     Tetraeder_storage* tet;
01547 
01548     sum=0.0;
01549     for(tet=bocedata->Give_tets();tet!=NULL;tet=tet->Next()) {
01550       if(tet->N0()==num_v) {
01551         sum = sum + ZM0*YM0 / (tet->Det()*6.0)
01552           * Coefficient_triangle;
01553       }
01554       else if(tet->N1()==num_v) {
01555         sum = sum + ZM1*YM1 / (tet->Det()*6.0)
01556           * Coefficient_triangle;
01557       }
01558       else if(tet->N2()==num_v) {
01559         sum = sum + ZM2*YM2 / (tet->Det()*6.0)
01560           * Coefficient_triangle;
01561       }
01562       else if(tet->N3()==num_v) {
01563         sum = sum + ZM3*YM3 / (tet->Det()*6.0)
01564           * Coefficient_triangle;
01565       }
01566     }     
01567     return sum;
01568   }
01569 };
01570 
01571 // Formeln fuer dzdy
01572 class dzdy_FE_variable {
01573  public:  
01574   diff_type typ_u() { return  Dz_type; }; 
01575   diff_type typ_v() { return  Dy_type; };
01576 
01577   dzdy_FE_variable() { }
01578   static inline int    ops_interior() { return 46; }  // 38 + 8
01579   static inline int    ops_diag_interior() { return 8; }
01580   static inline double stencil(double uN,
01581                                double uW, double uM, double uE,
01582                                double uS, double uT, double uD,
01583                                double uND, double uWN, double uWT,
01584                                double uED, double uST, double uES,
01585                                double uEST, double uWND,
01586                                double aN,
01587                                double aW, double aM, double aE,
01588                                double aS, double aT, double aD,
01589                                double aST, double aET, double aES,
01590                                double aND, double aWD, double aWN,
01591                                double aEST, double aWND,
01592                                double h_mesh) {
01593     return (
01594             var_interpo_WSD * (   uM - uD) +
01595             (var_interpo_ESD * (   uM - uD) +
01596              var_interpo_WNT * ( - uT + uM)) * 2.0 +
01597             var_interpo_WND * ( - uWN + uWND - uN + uND + uM - uD) +
01598             var_interpo_END * ( - uN + uND - uE + uED) +
01599             var_interpo_WST * (   uWT - uW + uST - uS) +
01600             var_interpo_EST * ( - uT + uM + uST - uS + uEST - uES) +
01601             var_interpo_ENT * ( - uT + uM)
01602             ) * h_mesh / 12.0; 
01603   }
01604   static inline double diag_stencil(double h_mesh,
01605                                     double aN,
01606                                     double aW, double aM, double aE,
01607                                     double aS, double aT, double aD,
01608                                     double aST, double aET, double aES,
01609                                     double aND, double aWD, double aWN,
01610                                     double aEST, double aWND) {
01611     return (
01612             var_interpo_WSD +
01613             (var_interpo_ESD +
01614              var_interpo_WNT) * 2.0 +
01615             var_interpo_WND +
01616             var_interpo_EST +
01617             var_interpo_ENT 
01618             ) * h_mesh / 12.0; 
01619   }
01620   static inline double Factor_reg_cell(double h_mesh) {
01621     return h_mesh / 6.0;
01622   }
01623   static inline double diag_F_reg_cell(dir_sons dir_v,
01624                                        double h_mesh, double* sten) {
01625     return sten[dir_v*9] * h_mesh / 6.0;
01626   }
01627 
01628   // Operator applied to variable with number
01629   static inline double F_bo_cell(const BoCeData* bocedata,
01630                                  int num_v, int num_variable,
01631                                  int num_coefficient) {
01632     static double sum;
01633     Tetraeder_storage* tet;
01634 
01635     sum=0.0;
01636     for(tet=bocedata->Give_tets();tet!=NULL;tet=tet->Next()) {
01637       if(tet->N0()==num_v) {
01638         sum = sum +
01639               YM0*(ZM0*UU0 + ZM1*UU1 + ZM2*UU2 + ZM3*UU3) / (tet->Det()*6.0)
01640           * Coefficient_triangle;
01641       }
01642       else if(tet->N1()==num_v) {
01643         sum = sum +
01644               YM1*(ZM0*UU0 + ZM1*UU1 + ZM2*UU2 + ZM3*UU3) / (tet->Det()*6.0)
01645           * Coefficient_triangle;
01646       }
01647       else if(tet->N2()==num_v) {
01648         sum = sum +
01649               YM2*(ZM0*UU0 + ZM1*UU1 + ZM2*UU2 + ZM3*UU3) / (tet->Det()*6.0)
01650           * Coefficient_triangle;
01651       }
01652       else if(tet->N3()==num_v) {
01653         sum = sum +
01654               YM3*(ZM0*UU0 + ZM1*UU1 + ZM2*UU2 + ZM3*UU3) / (tet->Det()*6.0)
01655           * Coefficient_triangle;
01656       }
01657     }     
01658     return sum;
01659   }
01660   // Operator applied to a vector U_array
01661   static inline double Local_matrix_bo_cell(const BoCeData* bocedata,
01662                                             double*  U_array, int num_v,
01663                                             int num_coefficient) {
01664     static double sum;
01665     Tetraeder_storage* tet;
01666 
01667     sum=0.0;
01668     for(tet=bocedata->Give_tets();tet!=NULL;tet=tet->Next()) {
01669       if(tet->N0()==num_v) {
01670         sum = sum +
01671           YM0*(ZM0*UAA0 + ZM1*UAA1 + ZM2*UAA2 + ZM3*UAA3) / (tet->Det()*6.0)
01672           * Coefficient_triangle;
01673       }
01674       else if(tet->N1()==num_v) {
01675         sum = sum +
01676           YM1*(ZM0*UAA0 + ZM1*UAA1 + ZM2*UAA2 + ZM3*UAA3) / (tet->Det()*6.0)
01677           * Coefficient_triangle;
01678       }
01679       else if(tet->N2()==num_v) {
01680         sum = sum +
01681           YM2*(ZM0*UAA0 + ZM1*UAA1 + ZM2*UAA2 + ZM3*UAA3) / (tet->Det()*6.0)
01682           * Coefficient_triangle;
01683       }
01684       else if(tet->N3()==num_v) {
01685         sum = sum +
01686           YM3*(ZM0*UAA0 + ZM1*UAA1 + ZM2*UAA2 + ZM3*UAA3) / (tet->Det()*6.0)
01687           * Coefficient_triangle;
01688       }
01689     }     
01690     return sum;
01691   }
01692   static inline double diag_F_bo_cell(const BoCeData* bocedata,
01693                                       int num_v, int num_coefficient) {
01694    static double sum;
01695     Tetraeder_storage* tet;
01696 
01697     sum=0.0;
01698     for(tet=bocedata->Give_tets();tet!=NULL;tet=tet->Next()) {
01699       if(tet->N0()==num_v) {
01700         sum = sum + YM0*ZM0 / (tet->Det()*6.0)
01701           * Coefficient_triangle;
01702       }
01703       else if(tet->N1()==num_v) {
01704         sum = sum + YM1*ZM1 / (tet->Det()*6.0)
01705           * Coefficient_triangle;
01706       }
01707       else if(tet->N2()==num_v) {
01708         sum = sum + YM2*ZM2 / (tet->Det()*6.0)
01709           * Coefficient_triangle;
01710       }
01711       else if(tet->N3()==num_v) {
01712         sum = sum + YM3*ZM3 / (tet->Det()*6.0)
01713           * Coefficient_triangle;
01714       }
01715     }     
01716     return sum;
01717   }
01718 };
01719 
01721 
01722 // Formeln fuer dxhelm
01723 class dxhelm_FE_variable {
01724  public:
01725   diff_type typ_u() { return  Dx_type; }; 
01726   diff_type typ_v() { return  Helm_type; };
01727 
01728   dxhelm_FE_variable() { }
01729   static inline int    ops_interior() { return 58; } // 50 + 8
01730   static inline int    ops_diag_interior() { return 19; }
01731   static inline double stencil(double uN,
01732                                double uW, double uM, double uE,
01733                                double uS, double uT, double uD,
01734                                double uND, double uWN, double uWT,
01735                                double uED, double uST, double uES,
01736                                double uEST, double uWND,
01737                                double aN,
01738                                double aW, double aM, double aE,
01739                                double aS, double aT, double aD,
01740                                double aST, double aET, double aES,
01741                                double aND, double aWD, double aWN,
01742                                double aEST, double aWND,
01743                                double h_mesh) {
01744     return (
01745             var_interpo_WSD * ( - uW + uM) +
01746             var_interpo_ESD * ( - uM - uD + uE + uED - uS + uES) +
01747             var_interpo_WND * ( - uWN + uN + 2.0*(uND - uWND - uW + uM)) +
01748             var_interpo_END * ( 2.0*(uE - uM) - uD + uED) +
01749             var_interpo_WST * ( uT - uWT + 2.0*(uM - uW)) +
01750             var_interpo_EST * ( 2.0*(uE - uM - uST + uEST) - uS + uES) +
01751             var_interpo_WNT * ( - uWN + uN - uWT - uW + uT + uM) +
01752             var_interpo_ENT * ( - uM + uE)
01753             ) * h_mesh * h_mesh / 48.0; 
01754   }
01755   static inline double diag_stencil(double h_mesh,
01756                                     double aN,
01757                                     double aW, double aM, double aE,
01758                                     double aS, double aT, double aD,
01759                                     double aST, double aET, double aES,
01760                                     double aND, double aWD, double aWN,
01761                                     double aEST, double aWND) {
01762     return (
01763             var_interpo_WSD -
01764             var_interpo_ESD +
01765             (var_interpo_WND -
01766              var_interpo_END +
01767              var_interpo_WST -
01768              var_interpo_EST) * 2.0 +
01769             var_interpo_WNT -
01770             var_interpo_ENT
01771             ) * h_mesh * h_mesh / 48.0; 
01772   }
01773   static inline double Factor_reg_cell(double h_mesh) {
01774     return h_mesh * h_mesh / 24.0;
01775   }
01776   static inline double diag_F_reg_cell(dir_sons dir_v,
01777                                        double h_mesh, double* sten) {
01778     return sten[dir_v*9] * h_mesh * h_mesh / 24.0;
01779   }
01780 
01781   // Operator applied to variable with number
01782   static inline double F_bo_cell(const BoCeData* bocedata,
01783                                  int num_v, int num_variable,
01784                                  int num_coefficient) {
01785     static double sum;
01786     Tetraeder_storage* tet;
01787 
01788     sum=0.0;
01789     for(tet=bocedata->Give_tets();tet!=NULL;tet=tet->Next()) {
01790       if(tet->N0()==num_v || tet->N1()==num_v  ||
01791          tet->N2()==num_v || tet->N3()==num_v ) {
01792         sum = sum +
01793               (XM0*UU0 + XM1*UU1 + XM2*UU2 + XM3*UU3) / 24.0 
01794           * Coefficient_triangle;
01795       }
01796     }
01797     return sum;
01798   }
01799   // Operator applied to a vector U_array
01800   static inline double Local_matrix_bo_cell(const BoCeData* bocedata,
01801                                             double*  U_array, int num_v,
01802                                             int num_coefficient) {
01803     static double sum;
01804     Tetraeder_storage* tet;
01805 
01806     sum=0.0;
01807     for(tet=bocedata->Give_tets();tet!=NULL;tet=tet->Next()) {
01808       if(tet->N0()==num_v || tet->N1()==num_v  ||
01809          tet->N2()==num_v || tet->N3()==num_v ) {
01810         sum = sum +
01811               (XM0*UAA0 + XM1*UAA1 + XM2*UAA2 + XM3*UAA3) / 24.0 
01812           * Coefficient_triangle;
01813       }
01814     }
01815     return sum;
01816   }
01817   static inline double diag_F_bo_cell(const BoCeData* bocedata,
01818                                       int num_v, int num_coefficient) {
01819    static double sum;
01820     Tetraeder_storage* tet;
01821 
01822     sum=0.0;
01823     for(tet=bocedata->Give_tets();tet!=NULL;tet=tet->Next()) {
01824       if(tet->N0()==num_v) sum = sum + XM0 / 24.0 
01825           * Coefficient_triangle;
01826       if(tet->N1()==num_v) sum = sum + XM1 / 24.0 
01827           * Coefficient_triangle;
01828       if(tet->N2()==num_v) sum = sum + XM2 / 24.0 
01829           * Coefficient_triangle;
01830       if(tet->N3()==num_v) sum = sum + XM3 / 24.0 
01831           * Coefficient_triangle;
01832     }
01833     return sum;
01834   }
01835 };
01836 
01837 // Formeln fuer helmdx
01838 class helmdx_FE_variable {
01839  public:
01840   diff_type typ_v() { return  Dx_type; }; 
01841   diff_type typ_u() { return  Helm_type; };
01842 
01843   helmdx_FE_variable() { }
01844   static inline int    ops_interior() { return 58; }
01845   static inline int    ops_diag_interior() { return 19; }
01846   static inline double stencil(double uN,
01847                                double uW, double uM, double uE,
01848                                double uS, double uT, double uD,
01849                                double uND, double uWN, double uWT,
01850                                double uED, double uST, double uES,
01851                                double uEST, double uWND,
01852                                double aN,
01853                                double aW, double aM, double aE,
01854                                double aS, double aT, double aD,
01855                                double aST, double aET, double aES,
01856                                double aND, double aWD, double aWN,
01857                                double aEST, double aWND,
01858                                double h_mesh) {
01859     return (
01860             var_interpo_WSD * (   uW + uM + uD + uS) +
01861             var_interpo_ESD * ( - uM - uE - uED - uES) +
01862             var_interpo_WND * (   uWN + 2.0*(uWND + uW + uM) + uD) +
01863             var_interpo_END * ( - uN - 2.0*(uND + uM + uE) - uED) +
01864             var_interpo_WST * (   uWT + 2.0*(uW + uM + uST) + uS) +
01865             var_interpo_EST * ( - uT - 2.0*(uM + uE + uEST) - uES) +
01866             var_interpo_WNT * (   uWN + uWT + uW + uM) +
01867             var_interpo_ENT * ( - uN - uT - uM - uE)
01868             ) * h_mesh * h_mesh / 48.0; 
01869   }
01870   static inline double diag_stencil(double h_mesh,
01871                                     double aN,
01872                                     double aW, double aM, double aE,
01873                                     double aS, double aT, double aD,
01874                                     double aST, double aET, double aES,
01875                                     double aND, double aWD, double aWN,
01876                                     double aEST, double aWND) {
01877     return (
01878             var_interpo_WSD -
01879             var_interpo_ESD +
01880             (var_interpo_WND -
01881              var_interpo_END +
01882              var_interpo_WST -
01883              var_interpo_EST) * 2.0 +
01884             var_interpo_WNT -
01885             var_interpo_ENT
01886             ) * h_mesh * h_mesh / 48.0; 
01887   }
01888   static inline double Factor_reg_cell(double h_mesh) {
01889     return h_mesh * h_mesh / 24.0;
01890   }
01891   static inline double diag_F_reg_cell(dir_sons dir_v,
01892                                        double h_mesh, double* sten) {
01893     return sten[dir_v*9] * h_mesh * h_mesh / 24.0;
01894   }
01895 
01896   // Operator applied to variable with number
01897   static inline double F_bo_cell(const BoCeData* bocedata,
01898                                  int num_v, int num_variable,
01899                                  int num_coefficient) {
01900     static double sum;
01901     Tetraeder_storage* tet;
01902 
01903     sum=0.0;
01904     
01905     for(tet=bocedata->Give_tets();tet!=NULL;tet=tet->Next()) {
01906       if(tet->N0()==num_v) sum += XM0*(UU0 + UU1 + UU2 + UU3) / 24.0 
01907           * Coefficient_triangle;
01908       if(tet->N1()==num_v) sum += XM1*(UU0 + UU1 + UU2 + UU3) / 24.0 
01909           * Coefficient_triangle;
01910       if(tet->N2()==num_v) sum += XM2*(UU0 + UU1 + UU2 + UU3) / 24.0 
01911           * Coefficient_triangle;
01912       if(tet->N3()==num_v) sum += XM3*(UU0 + UU1 + UU2 + UU3) / 24.0 
01913           * Coefficient_triangle;
01914     }
01915     return sum;
01916   }
01917   // Operator applied to a vector U_array
01918   static inline double Local_matrix_bo_cell(const BoCeData* bocedata,
01919                                             double*  U_array, int num_v,
01920                                             int num_coefficient) {
01921    static double sum;
01922     Tetraeder_storage* tet;
01923 
01924     sum=0.0;
01925     
01926     for(tet=bocedata->Give_tets();tet!=NULL;tet=tet->Next()) {
01927       if(tet->N0()==num_v) sum += XM0*(UAA0 + UAA1 + UAA2 + UAA3) / 24.0 
01928           * Coefficient_triangle;
01929       if(tet->N1()==num_v) sum += XM1*(UAA0 + UAA1 + UAA2 + UAA3) / 24.0 
01930           * Coefficient_triangle;
01931       if(tet->N2()==num_v) sum += XM2*(UAA0 + UAA1 + UAA2 + UAA3) / 24.0 
01932           * Coefficient_triangle;
01933       if(tet->N3()==num_v) sum += XM3*(UAA0 + UAA1 + UAA2 + UAA3) / 24.0 
01934           * Coefficient_triangle;
01935     }
01936     return sum;
01937   }
01938   static inline double diag_F_bo_cell(const BoCeData* bocedata,
01939                                       int num_v, int num_coefficient) {
01940    static double sum;
01941     Tetraeder_storage* tet;
01942 
01943     sum=0.0;
01944     
01945     for(tet=bocedata->Give_tets();tet!=NULL;tet=tet->Next()) {
01946       if(tet->N0()==num_v) sum += XM0 / 24.0 
01947           * Coefficient_triangle;
01948       if(tet->N1()==num_v) sum += XM1 / 24.0 
01949           * Coefficient_triangle;
01950       if(tet->N2()==num_v) sum += XM2 / 24.0 
01951           * Coefficient_triangle;
01952       if(tet->N3()==num_v) sum += XM3 / 24.0 
01953           * Coefficient_triangle;
01954     }
01955     return sum;
01956   }
01957 };
01958 
01959 // Formeln fuer dyhelm
01960 class dyhelm_FE_variable {
01961  public:
01962   diff_type typ_u() { return  Dy_type; }; 
01963   diff_type typ_v() { return  Helm_type; };
01964 
01965   dyhelm_FE_variable() { }
01966   static inline int    ops_interior() { return 68; }
01967   static inline int    ops_diag_interior() { return 19; }
01968   static inline double stencil(double uN,
01969                                double uW, double uM, double uE,
01970                                double uS, double uT, double uD,
01971                                double uND, double uWN, double uWT,
01972                                double uED, double uST, double uES,
01973                                double uEST, double uWND,
01974                                double aN,
01975                                double aW, double aM, double aE,
01976                                double aS, double aT, double aD,
01977                                double aST, double aET, double aES,
01978                                double aND, double aWD, double aWN,
01979                                double aEST, double aWND,
01980                                double h_mesh) {
01981     return (
01982             var_interpo_WSD * (   uM - uS) +
01983             var_interpo_ESD * (   2.0*uM - uD + uE + uED - uS - 2.0*uES) +
01984             var_interpo_WND * ( uN - uM + 2.0*(uWN + uND - uW - uD)) +
01985             var_interpo_END * (   uN + 2.0*uND - 2.0*uM - uD + uE - uED) +
01986             var_interpo_WST * (   uWT - uW + uT + 2.0*uM - 2.0*uST - uS) +
01987             var_interpo_EST * ( uM - uS + 2.0*(uT + uE - uST - uES)) +
01988             var_interpo_WNT * (   2.0*uWN + uN - uWT - uW + uT - 2.0*uM) +
01989             var_interpo_ENT * (   uN - uM)
01990             ) * h_mesh * h_mesh / 48.0; 
01991   }
01992   static inline double diag_stencil(double h_mesh,
01993                                     double aN,
01994                                     double aW, double aM, double aE,
01995                                     double aS, double aT, double aD,
01996                                     double aST, double aET, double aES,
01997                                     double aND, double aWD, double aWN,
01998                                     double aEST, double aWND) {
01999     return (
02000             var_interpo_WSD +
02001             (var_interpo_ESD -
02002             var_interpo_END +
02003             var_interpo_WST -
02004             var_interpo_WNT) * 2.0 -
02005             var_interpo_WND +
02006             var_interpo_EST -
02007             var_interpo_ENT 
02008             ) * h_mesh * h_mesh / 48.0; 
02009   }
02010   static inline double Factor_reg_cell(double h_mesh) {
02011     return h_mesh * h_mesh / 24.0;
02012   }
02013   static inline double diag_F_reg_cell(dir_sons dir_v,
02014                                        double h_mesh, double* sten) {
02015     return sten[dir_v*9] * h_mesh * h_mesh / 24.0;
02016   }
02017 
02018   // Operator applied to variable with number
02019   static inline double F_bo_cell(const BoCeData* bocedata,
02020                                  int num_v, int num_variable,
02021                                  int num_coefficient) {
02022     static double sum;
02023     Tetraeder_storage* tet;
02024 
02025     sum=0.0;
02026     for(tet=bocedata->Give_tets();tet!=NULL;tet=tet->Next()) {
02027       if(tet->N0()==num_v || tet->N1()==num_v  ||
02028          tet->N2()==num_v || tet->N3()==num_v ) {
02029         sum = sum +
02030               (YM0*UU0 + YM1*UU1 + YM2*UU2 + YM3*UU3) / 24.0 
02031           * Coefficient_triangle;
02032       }
02033     }
02034     //    return 0.0;
02035     return sum;
02036   }
02037   // Operator applied to a vector U_array
02038   static inline double Local_matrix_bo_cell(const BoCeData* bocedata,
02039                                             double*  U_array, int num_v,
02040                                             int num_coefficient) {
02041    static double sum;
02042     Tetraeder_storage* tet;
02043 
02044     sum=0.0;
02045     for(tet=bocedata->Give_tets();tet!=NULL;tet=tet->Next()) {
02046       if(tet->N0()==num_v || tet->N1()==num_v  ||
02047          tet->N2()==num_v || tet->N3()==num_v ) {
02048         sum = sum +
02049               (YM0*UAA0 + YM1*UAA1 + YM2*UAA2 + YM3*UAA3) / 24.0 
02050           * Coefficient_triangle;
02051       }
02052     }
02053     //    return 0.0;
02054     return sum;
02055   }
02056   static inline double diag_F_bo_cell(const BoCeData* bocedata,
02057                                       int num_v, int num_coefficient) {
02058    static double sum;
02059     Tetraeder_storage* tet;
02060 
02061     sum=0.0;
02062     for(tet=bocedata->Give_tets();tet!=NULL;tet=tet->Next()) {
02063       if(tet->N0()==num_v) sum = sum + YM0 / 24.0 
02064           * Coefficient_triangle;
02065       if(tet->N1()==num_v) sum = sum + YM1 / 24.0 
02066           * Coefficient_triangle;
02067       if(tet->N2()==num_v) sum = sum + YM2 / 24.0 
02068           * Coefficient_triangle;
02069       if(tet->N3()==num_v) sum = sum + YM3 / 24.0 
02070           * Coefficient_triangle;
02071     }
02072     return sum;
02073   }
02074 };
02075 
02076 // Formeln fuer helmdy
02077 class helmdy_FE_variable {
02078  public:
02079   diff_type typ_v() { return  Dy_type; }; 
02080   diff_type typ_u() { return  Helm_type; };
02081 
02082   helmdy_FE_variable() { }
02083   static inline int    ops_interior() { return 54; } // 46 + 8
02084   static inline int    ops_diag_interior() { return 19; }
02085   static inline double stencil(double uN,
02086                                double uW, double uM, double uE,
02087                                double uS, double uT, double uD,
02088                                double uND, double uWN, double uWT,
02089                                double uED, double uST, double uES,
02090                                double uEST, double uWND,
02091                                double aN,
02092                                double aW, double aM, double aE,
02093                                double aS, double aT, double aD,
02094                                double aST, double aET, double aES,
02095                                double aND, double aWD, double aWN,
02096                                double aEST, double aWND,
02097                                double h_mesh) {
02098     return (
02099             var_interpo_WSD * (   uW + uM + uD + uS) +
02100             var_interpo_ESD * (   2.0*(uM + uD + uES) + uED + uS) +
02101             var_interpo_WND * ( - 2.0*(uWN + uND) - uN + uW - uM + uD) +
02102             var_interpo_END * ( - uN - 2.0*(uND + uM + uE) - uED) +
02103             var_interpo_WST * (   uWT + 2.0*(uW + uM + uST) + uS) +
02104             var_interpo_EST * ( - uT + uM - uE + 2.0*(uST + uES) + uS) +
02105             var_interpo_WNT * ( - 2.0*(uWN + uT + uM) - uN - uWT) +
02106             var_interpo_ENT * ( - uN - uT - uM - uE)
02107             ) * h_mesh * h_mesh / 48.0; 
02108   }
02109   static inline double diag_stencil(double h_mesh,
02110                                     double aN,
02111                                     double aW, double aM, double aE,
02112                                     double aS, double aT, double aD,
02113                                     double aST, double aET, double aES,
02114                                     double aND, double aWD, double aWN,
02115                                     double aEST, double aWND) {
02116     return (
02117             var_interpo_WSD +
02118             (var_interpo_ESD -
02119              var_interpo_END -
02120              var_interpo_WNT +
02121              var_interpo_WST) * 2.0 -
02122             var_interpo_WND +
02123             var_interpo_EST -
02124             var_interpo_ENT 
02125             ) * h_mesh * h_mesh / 48.0; 
02126   }
02127   static inline double Factor_reg_cell(double h_mesh) {
02128     return h_mesh * h_mesh / 24.0;
02129   }
02130   static inline double diag_F_reg_cell(dir_sons dir_v,
02131                                        double h_mesh, double* sten) {
02132     return sten[dir_v*9] * h_mesh * h_mesh / 24.0; 
02133   }
02134 
02135   // Operator applied to variable with number
02136   static inline double F_bo_cell(const BoCeData* bocedata,
02137                                  int num_v, int num_variable,
02138                                  int num_coefficient) {
02139     static double sum;
02140     Tetraeder_storage* tet;
02141 
02142     sum=0.0;
02143     for(tet=bocedata->Give_tets();tet!=NULL;tet=tet->Next()) {
02144       if(tet->N0()==num_v) sum += YM0*(UU0 + UU1 + UU2 + UU3) / 24.0 
02145           * Coefficient_triangle;
02146       if(tet->N1()==num_v) sum += YM1*(UU0 + UU1 + UU2 + UU3) / 24.0 
02147           * Coefficient_triangle;
02148       if(tet->N2()==num_v) sum += YM2*(UU0 + UU1 + UU2 + UU3) / 24.0 
02149           * Coefficient_triangle;
02150       if(tet->N3()==num_v) sum += YM3*(UU0 + UU1 + UU2 + UU3) / 24.0 
02151           * Coefficient_triangle;
02152     }
02153     return sum;
02154   }
02155   // Operator applied to a vector U_array
02156   static inline double Local_matrix_bo_cell(const BoCeData* bocedata,
02157                                             double*  U_array, int num_v,
02158                                             int num_coefficient) {
02159    static double sum;
02160     Tetraeder_storage* tet;
02161 
02162     sum=0.0;
02163     for(tet=bocedata->Give_tets();tet!=NULL;tet=tet->Next()) {
02164       if(tet->N0()==num_v) sum += YM0*(UAA0 + UAA1 + UAA2 + UAA3) / 24.0 
02165           * Coefficient_triangle;
02166       if(tet->N1()==num_v) sum += YM1*(UAA0 + UAA1 + UAA2 + UAA3) / 24.0 
02167           * Coefficient_triangle;
02168       if(tet->N2()==num_v) sum += YM2*(UAA0 + UAA1 + UAA2 + UAA3) / 24.0 
02169           * Coefficient_triangle;
02170       if(tet->N3()==num_v) sum += YM3*(UAA0 + UAA1 + UAA2 + UAA3) / 24.0 
02171           * Coefficient_triangle;
02172     }
02173     return sum;
02174   }
02175   static inline double diag_F_bo_cell(const BoCeData* bocedata,
02176                                       int num_v, int num_coefficient) {
02177    static double sum;
02178     Tetraeder_storage* tet;
02179 
02180     sum=0.0;
02181     for(tet=bocedata->Give_tets();tet!=NULL;tet=tet->Next()) {
02182       if(tet->N0()==num_v) sum += YM0 / 24.0 
02183           * Coefficient_triangle;
02184       if(tet->N1()==num_v) sum += YM1 / 24.0 
02185           * Coefficient_triangle;
02186       if(tet->N2()==num_v) sum += YM2 / 24.0 
02187           * Coefficient_triangle;
02188       if(tet->N3()==num_v) sum += YM3 / 24.0 
02189           * Coefficient_triangle;
02190     }
02191     return sum;
02192   }
02193 };
02194 
02195 // Formeln fuer dzhelm
02196 class dzhelm_FE_variable {
02197  public:
02198   diff_type typ_u() { return  Dz_type; }; 
02199   diff_type typ_v() { return  Helm_type; };
02200 
02201   dzhelm_FE_variable() { }
02202   static inline int    ops_interior() { return 58; }
02203   static inline int    ops_diag_interior() { return 19; }
02204   static inline double stencil(double uN,
02205                                double uW, double uM, double uE,
02206                                double uS, double uT, double uD,
02207                                double uND, double uWN, double uWT,
02208                                double uED, double uST, double uES,
02209                                double uEST, double uWND,
02210                                double aN,
02211                                double aW, double aM, double aE,
02212                                double aS, double aT, double aD,
02213                                double aST, double aET, double aES,
02214                                double aND, double aWD, double aWN,
02215                                double aEST, double aWND,
02216                                double h_mesh) {
02217     return (
02218             var_interpo_WSD * (   uM - uD) +
02219             var_interpo_ESD * (   2.0*(uM - uD) + uE - uED) +
02220             var_interpo_WND * ( uN - uND + 2.0*(uWN - uWND + uM - uD)) +
02221             var_interpo_END * (   uN - uND + uM - uD + uE - uED) +
02222             var_interpo_WST * (   uWT - uW + uT - uM + uST - uS) +
02223             var_interpo_EST * ( uST - uS + 2.0*(uT - uM + uEST - uES)) +
02224             var_interpo_WNT * (   uWT - uW + 2.0*(uT - uM)) +
02225             var_interpo_ENT * (   uT - uM) 
02226             )  * h_mesh * h_mesh / 48.0; 
02227   }
02228   static inline double diag_stencil(double h_mesh,
02229                                     double aN,
02230                                     double aW, double aM, double aE,
02231                                     double aS, double aT, double aD,
02232                                     double aST, double aET, double aES,
02233                                     double aND, double aWD, double aWN,
02234                                     double aEST, double aWND) {
02235     return (
02236             var_interpo_WSD +
02237             (var_interpo_ESD -
02238              var_interpo_EST -
02239              var_interpo_WNT +
02240              var_interpo_WND) * 2.0 +
02241             var_interpo_END -
02242             var_interpo_WST -
02243             var_interpo_ENT
02244             )  * h_mesh * h_mesh / 48.0; 
02245   }
02246   static inline double Factor_reg_cell(double h_mesh) {
02247     return h_mesh * h_mesh / 24.0;
02248   }
02249   static inline double diag_F_reg_cell(dir_sons dir_v,
02250                                        double h_mesh, double* sten) {
02251     return sten[dir_v*9] * h_mesh * h_mesh / 24.0;
02252   }
02253 
02254   // Operator applied to variable with number
02255   static inline double F_bo_cell(const BoCeData* bocedata,
02256                                  int num_v, int num_variable,
02257                                  int num_coefficient) {
02258     static double sum;
02259     Tetraeder_storage* tet;
02260 
02261     sum=0.0;
02262     for(tet=bocedata->Give_tets();tet!=NULL;tet=tet->Next()) {
02263       if(tet->N0()==num_v || tet->N1()==num_v  ||
02264          tet->N2()==num_v || tet->N3()==num_v ) {
02265         sum = sum +
02266               (ZM0*UU0 + ZM1*UU1 + ZM2*UU2 + ZM3*UU3) / 24.0 
02267           * Coefficient_triangle;
02268       }
02269     }
02270     // return 0.0;
02271     return sum;
02272   }
02273   // Operator applied to a vector U_array
02274   static inline double Local_matrix_bo_cell(const BoCeData* bocedata,
02275                                             double*  U_array, int num_v,
02276                                             int num_coefficient) {
02277   static double sum;
02278     Tetraeder_storage* tet;
02279 
02280     sum=0.0;
02281     for(tet=bocedata->Give_tets();tet!=NULL;tet=tet->Next()) {
02282       if(tet->N0()==num_v || tet->N1()==num_v  ||
02283          tet->N2()==num_v || tet->N3()==num_v ) {
02284         sum = sum +
02285               (ZM0*UAA0 + ZM1*UAA1 + ZM2*UAA2 + ZM3*UAA3) / 24.0 
02286           * Coefficient_triangle;
02287       }
02288     }
02289     // return 0.0;
02290     return sum;
02291   }
02292   static inline double diag_F_bo_cell(const BoCeData* bocedata,
02293                                       int num_v, int num_coefficient) {
02294    static double sum;
02295     Tetraeder_storage* tet;
02296 
02297     sum=0.0;
02298     for(tet=bocedata->Give_tets();tet!=NULL;tet=tet->Next()) {
02299       if(tet->N0()==num_v) sum = sum + ZM0 / 24.0 
02300           * Coefficient_triangle;
02301       if(tet->N1()==num_v) sum = sum + ZM1 / 24.0 
02302           * Coefficient_triangle;
02303       if(tet->N2()==num_v) sum = sum + ZM2 / 24.0 
02304           * Coefficient_triangle;
02305       if(tet->N3()==num_v) sum = sum + ZM3 / 24.0 
02306           * Coefficient_triangle;
02307     }
02308     return sum;
02309   }
02310 };
02311 
02312 // Formeln fuer helmdz
02313 class helmdz_FE_variable {
02314  public:
02315   diff_type typ_v() { return  Dz_type; }; 
02316   diff_type typ_u() { return  Helm_type; };
02317 
02318   helmdz_FE_variable() { }
02319   static inline int    ops_interior() { return 58; }
02320   static inline int    ops_diag_interior() { return 19; }
02321   static inline double stencil(double uN,
02322                                double uW, double uM, double uE,
02323                                double uS, double uT, double uD,
02324                                double uND, double uWN, double uWT,
02325                                double uED, double uST, double uES,
02326                                double uEST, double uWND,
02327                                double aN,
02328                                double aW, double aM, double aE,
02329                                double aS, double aT, double aD,
02330                                double aST, double aET, double aES,
02331                                double aND, double aWD, double aWN,
02332                                double aEST, double aWND,
02333                                double h_mesh) {
02334     return (
02335             var_interpo_WSD * (   uW + uM + uD + uS) +
02336             var_interpo_ESD * (   2.0*(uM + uD + uES) + uED + uS) +
02337             var_interpo_WND * (   uND + uW + 2.0*(uWND + uM + uD)) +
02338             var_interpo_END * (   uND + uM + uD + uED) +
02339             var_interpo_WST * ( - uWT - uT - uM - uST) +
02340             var_interpo_EST * ( - uE - uST - 2.0*(uT + uM + uEST)) +
02341             var_interpo_WNT * ( - uN - uWT - 2.0*(uWN + uT + uM)) +
02342             var_interpo_ENT * ( - uN - uT - uM - uE) 
02343             ) * h_mesh * h_mesh / 48.0; 
02344   }
02345   static inline double diag_stencil(double h_mesh,
02346                                     double aN,
02347                                     double aW, double aM, double aE,
02348                                     double aS, double aT, double aD,
02349                                     double aST, double aET, double aES,
02350                                     double aND, double aWD, double aWN,
02351                                     double aEST, double aWND) {
02352     return (
02353             var_interpo_WSD +
02354             (var_interpo_ESD -
02355              var_interpo_EST -
02356              var_interpo_WNT +
02357              var_interpo_WND) * 2.0 +
02358             var_interpo_END -
02359             var_interpo_WST -
02360             var_interpo_ENT
02361             )  * h_mesh * h_mesh / 48.0; 
02362   }
02363   static inline double Factor_reg_cell(double h_mesh) {
02364     return h_mesh * h_mesh / 24.0;
02365   }
02366   static inline double diag_F_reg_cell(dir_sons dir_v,
02367                                        double h_mesh, double* sten) {
02368     return sten[dir_v*9] * h_mesh * h_mesh / 24.0;
02369   }
02370 
02371   // Operator applied to variable with number
02372   static inline double F_bo_cell(const BoCeData* bocedata,
02373                                  int num_v, int num_variable,
02374                                  int num_coefficient) {
02375     static double sum;
02376     Tetraeder_storage* tet;
02377 
02378     sum=0.0;
02379     for(tet=bocedata->Give_tets();tet!=NULL;tet=tet->Next()) {
02380       if(tet->N0()==num_v) sum += ZM0*(UU0 + UU1 + UU2 + UU3) / 24.0 
02381           * Coefficient_triangle;
02382       if(tet->N1()==num_v) sum += ZM1*(UU0 + UU1 + UU2 + UU3) / 24.0 
02383           * Coefficient_triangle;
02384       if(tet->N2()==num_v) sum += ZM2*(UU0 + UU1 + UU2 + UU3) / 24.0 
02385           * Coefficient_triangle;
02386       if(tet->N3()==num_v) sum += ZM3*(UU0 + UU1 + UU2 + UU3) / 24.0 
02387           * Coefficient_triangle;
02388     }
02389     // return 0.0;
02390     return sum;
02391   }
02392   // Operator applied to a vector U_array
02393   static inline double Local_matrix_bo_cell(const BoCeData* bocedata,
02394                                             double*  U_array, int num_v,
02395                                             int num_coefficient) {
02396    static double sum;
02397     Tetraeder_storage* tet;
02398 
02399     sum=0.0;
02400     for(tet=bocedata->Give_tets();tet!=NULL;tet=tet->Next()) {
02401       if(tet->N0()==num_v) sum += ZM0*(UAA0 + UAA1 + UAA2 + UAA3) / 24.0 
02402           * Coefficient_triangle;
02403       if(tet->N1()==num_v) sum += ZM1*(UAA0 + UAA1 + UAA2 + UAA3) / 24.0 
02404           * Coefficient_triangle;
02405       if(tet->N2()==num_v) sum += ZM2*(UAA0 + UAA1 + UAA2 + UAA3) / 24.0 
02406           * Coefficient_triangle;
02407       if(tet->N3()==num_v) sum += ZM3*(UAA0 + UAA1 + UAA2 + UAA3) / 24.0 
02408           * Coefficient_triangle;
02409     }
02410     // return 0.0;
02411     return sum;
02412   }
02413   static inline double diag_F_bo_cell(const BoCeData* bocedata,
02414                                       int num_v, int num_coefficient) {
02415     static double sum;
02416     Tetraeder_storage* tet;
02417 
02418     sum=0.0;
02419     for(tet=bocedata->Give_tets();tet!=NULL;tet=tet->Next()) {
02420       if(tet->N0()==num_v) sum += ZM0 / 24.0 
02421           * Coefficient_triangle;
02422       if(tet->N1()==num_v) sum += ZM1 / 24.0 
02423           * Coefficient_triangle;
02424       if(tet->N2()==num_v) sum += ZM2 / 24.0 
02425           * Coefficient_triangle;
02426       if(tet->N3()==num_v) sum += ZM3 / 24.0 
02427           * Coefficient_triangle;
02428     }
02429     return sum;
02430   }
02431 };
02432 
02433 
02434 #endif
02435 
02436