Main Page | Namespace List | Class Hierarchy | Class List | File List | Class Members | File Members

src/expde/formulas/diffopcv.h

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