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

src/expde/formulas/diffopc.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 // diffopc.h
00019 //
00020 // ------------------------------------------------------------
00021 
00022 #ifndef DIFFOPC_H_
00023 #define DIFFOPC_H_
00024        
00025 // -----------------------------------------------------------
00026 //    FFFF    OO     RR     M     M   EEEE  L     N     N
00027 //    F      O  O    R R    M M M M   E     L     N N   N
00028 //    FFF    O  O    RR     M  M  M   EEEE  L     N  N  N
00029 //    F      O  O    R R    M     M   E     L     N   N N
00030 //    F       OO     R  R   M     M   EEEE  LLLL  N     N
00031 // -----------------------------------------------------------
00032 
00034 // 3.:  Constant coefficient Differential operators
00036 
00037 // Formeln fuer den Laplace Operator
00038 
00039 class laplace_FE_const {
00040  public:
00041   diff_type typ_u() { return  Poisson_type; }; 
00042   diff_type typ_v() { return  Poisson_type; };
00043 
00044   laplace_FE_const() { }
00045   static inline int    ops_interior() { return 19; }
00046   static inline int    ops_diag_interior() { return 1; }
00047   static inline double stencil(double uN,
00048                                double uW, double uM, double uO,
00049                                double uS, double uT, double uD,
00050                                double uND, double uWN, double uWT,
00051                                double uED, double uST, double uES,
00052                                double uEST, double uWND,
00053                                double h_mesh) {
00054     return  (20.0 * uM - 4.0*(uT + uD + uO + uW) - 2.0*(uN + uS)
00055                  -  (uND + uWN - uWT - uED + uST + uES -
00056                     uEST - uWND)) * h_mesh / 3.0; 
00057   }
00058   static inline double diag_stencil(double h_mesh) {
00059     return 20.0 / 3.0 * h_mesh;
00060   }
00061   static inline double F_reg_cell(double*const * u_Recell, int num_var,
00062                                   dir_sons dir_v,
00063                                   double h_mesh, double* sten) {
00064     static int i;
00065     static double sum;
00066     sum = 0.0;
00067     for(i=0;i<8;++i) sum = sum + sten[dir_v+8*i] * u_Recell[i][num_var];
00068     return sum * h_mesh / 6.0;
00069   }
00070   static inline double Factor_reg_cell(double h_mesh) {
00071     return h_mesh / 6.0;
00072   }
00073   static inline double diag_F_reg_cell(dir_sons dir_v,
00074                                        double h_mesh, double* sten) {
00075     return sten[dir_v*9] * h_mesh / 6.0;
00076   }
00077 
00078   // Operator applied to variable with number: num_variable
00079   // with test function: num_v
00080   // on a boundary cell
00081   static inline double F_bo_cell(const BoCeData* bocedata,
00082                                  int num_v, int num_variable) {
00083     static double sum;
00084     Tetraeder_storage* tet;
00085     
00086 
00087     sum=0.0;
00088     for(tet=bocedata->Give_tets();tet!=NULL;tet=tet->Next()) {
00089       if(tet->N0()==num_v) {
00090         sum = sum + (
00091               XM0*(XM0*UU0 + XM1*UU1 + XM2*UU2 + XM3*UU3) +
00092               YM0*(YM0*UU0 + YM1*UU1 + YM2*UU2 + YM3*UU3) +
00093               ZM0*(ZM0*UU0 + ZM1*UU1 + ZM2*UU2 + ZM3*UU3)) / (tet->Det()*6.0);
00094       }
00095       else if(tet->N1()==num_v) {
00096         sum = sum + (
00097               XM1*(XM0*UU0 + XM1*UU1 + XM2*UU2 + XM3*UU3) +
00098               YM1*(YM0*UU0 + YM1*UU1 + YM2*UU2 + YM3*UU3) +
00099               ZM1*(ZM0*UU0 + ZM1*UU1 + ZM2*UU2 + ZM3*UU3)) / (tet->Det()*6.0);
00100       }
00101       else if(tet->N2()==num_v) {
00102         sum = sum + (
00103               XM2*(XM0*UU0 + XM1*UU1 + XM2*UU2 + XM3*UU3) +
00104               YM2*(YM0*UU0 + YM1*UU1 + YM2*UU2 + YM3*UU3) +
00105               ZM2*(ZM0*UU0 + ZM1*UU1 + ZM2*UU2 + ZM3*UU3)) / (tet->Det()*6.0);
00106       }
00107       else if(tet->N3()==num_v) {
00108         sum = sum + (
00109               XM3*(XM0*UU0 + XM1*UU1 + XM2*UU2 + XM3*UU3) +
00110               YM3*(YM0*UU0 + YM1*UU1 + YM2*UU2 + YM3*UU3) +
00111               ZM3*(ZM0*UU0 + ZM1*UU1 + ZM2*UU2 + ZM3*UU3)) / (tet->Det()*6.0);
00112       }
00113     }
00114     return sum;
00115   }
00116 
00117   // Operator applied to a vector U_array
00118   // with test function: num_v
00119   // on a boundary cell
00120   static inline double Local_matrix_bo_cell(const BoCeData* bocedata,
00121                                           double*  U_array, int num_v) {
00122     static double sum;
00123     Tetraeder_storage* tet;
00124     
00125     sum=0.0;
00126     for(tet=bocedata->Give_tets();tet!=NULL;tet=tet->Next()) {
00127       if(tet->N0()==num_v) {
00128         sum = sum + (
00129                      XM0*(XM0*UAA0 + XM1*UAA1 + XM2*UAA2 + XM3*UAA3) +
00130                      YM0*(YM0*UAA0 + YM1*UAA1 + YM2*UAA2 + YM3*UAA3) +
00131                      ZM0*(ZM0*UAA0 + ZM1*UAA1 + ZM2*UAA2 + ZM3*UAA3)) /
00132           (tet->Det()*6.0);
00133       }
00134       else if(tet->N1()==num_v) {
00135         sum = sum + (
00136                      XM1*(XM0*UAA0 + XM1*UAA1 + XM2*UAA2 + XM3*UAA3) +
00137                      YM1*(YM0*UAA0 + YM1*UAA1 + YM2*UAA2 + YM3*UAA3) +
00138                      ZM1*(ZM0*UAA0 + ZM1*UAA1 + ZM2*UAA2 + ZM3*UAA3)) /
00139           (tet->Det()*6.0);
00140       }
00141       else if(tet->N2()==num_v) {
00142         sum = sum + (
00143                      XM2*(XM0*UAA0 + XM1*UAA1 + XM2*UAA2 + XM3*UAA3) +
00144                      YM2*(YM0*UAA0 + YM1*UAA1 + YM2*UAA2 + YM3*UAA3) +
00145                      ZM2*(ZM0*UAA0 + ZM1*UAA1 + ZM2*UAA2 + ZM3*UAA3)) /
00146           (tet->Det()*6.0);
00147       }
00148       else if(tet->N3()==num_v) {
00149         sum = sum + (
00150                      XM3*(XM0*UAA0 + XM1*UAA1 + XM2*UAA2 + XM3*UAA3) +
00151                      YM3*(YM0*UAA0 + YM1*UAA1 + YM2*UAA2 + YM3*UAA3) +
00152                      ZM3*(ZM0*UAA0 + ZM1*UAA1 + ZM2*UAA2 + ZM3*UAA3)) /
00153           (tet->Det()*6.0);
00154       }
00155     }
00156     return sum;  
00157   }
00158 
00159   static inline double diag_F_bo_cell(const BoCeData* bocedata,
00160                                       int num_v) {
00161     static double sum;
00162     Tetraeder_storage* tet;
00163     
00164     sum=0.0;
00165     
00166     for(tet=bocedata->Give_tets();tet!=NULL;tet=tet->Next()) {
00167       if(tet->N0()==num_v) {
00168         sum = sum + (XM0*XM0 + YM0*YM0 + ZM0*ZM0) / (tet->Det()*6.0);
00169       }
00170       else if(tet->N1()==num_v) {
00171         sum = sum + (XM1*XM1 + YM1*YM1 + ZM1*ZM1) / (tet->Det()*6.0);
00172       }
00173       else if(tet->N2()==num_v) {
00174         sum = sum + (XM2*XM2 + YM2*YM2 + ZM2*ZM2) / (tet->Det()*6.0);
00175       }
00176       else if(tet->N3()==num_v) {
00177         sum = sum + (XM3*XM3 + YM3*YM3 + ZM3*ZM3) / (tet->Det()*6.0);
00178       }
00179     }
00180     return sum;
00181   }
00182 };
00183 
00184 
00185 
00186 // Formeln fuer dxdx
00187 class dxdx_FE_const {
00188  public:
00189   diff_type typ_u() { return  Dx_type; }; 
00190   diff_type typ_v() { return  Dx_type; };
00191 
00192   dxdx_FE_const() { }
00193   static inline int    ops_interior() { return 4; }
00194   static inline int    ops_diag_interior() { return 1; }
00195   static inline double stencil(double uN,
00196                                double uW, double uM, double uO,
00197                                double uS, double uT, double uD,
00198                                double uND, double uWN, double uWT,
00199                                double uED, double uST, double uES,
00200                                double uEST, double uWND,
00201                                double h_mesh) {
00202     return  (  2.0 * uM - uW - uO) * h_mesh; 
00203   }
00204   static inline double diag_stencil(double h_mesh) {
00205     return 2.0 * h_mesh;
00206   }
00207   static inline double F_reg_cell(double*const * u_Recell, int num_var,
00208                                   dir_sons dir_v,
00209                                   double h_mesh, double* sten) {
00210     static int i;
00211     static double sum;
00212   
00213     sum = 0.0;
00214     for(i=0;i<8;++i) sum = sum + sten[dir_v+8*i] * u_Recell[i][num_var];
00215     return sum * h_mesh / 6.0;
00216   }
00217   static inline double Factor_reg_cell(double h_mesh) {
00218     return h_mesh / 6.0;
00219   }
00220   static inline double diag_F_reg_cell(dir_sons dir_v,
00221                                        double h_mesh, double* sten) {
00222     return sten[dir_v*9] * h_mesh / 6.0;
00223   }
00224 
00225   // Operator applied to variable with number
00226   static inline double F_bo_cell(const BoCeData* bocedata,
00227                                  int num_v, int num_variable) {
00228     static double sum;
00229     Tetraeder_storage* tet;
00230 
00231     sum=0.0;
00232     for(tet=bocedata->Give_tets();tet!=NULL;tet=tet->Next()) {
00233       if(tet->N0()==num_v) {
00234         sum = sum +
00235               XM0*(XM0*UU0 + XM1*UU1 + XM2*UU2 + XM3*UU3) / (tet->Det()*6.0);
00236       }
00237       else if(tet->N1()==num_v) {
00238         sum = sum +
00239               XM1*(XM0*UU0 + XM1*UU1 + XM2*UU2 + XM3*UU3) / (tet->Det()*6.0);
00240       }
00241       else if(tet->N2()==num_v) {
00242         sum = sum +
00243               XM2*(XM0*UU0 + XM1*UU1 + XM2*UU2 + XM3*UU3) / (tet->Det()*6.0);
00244       }
00245       else if(tet->N3()==num_v) {
00246         sum = sum +
00247               XM3*(XM0*UU0 + XM1*UU1 + XM2*UU2 + XM3*UU3) / (tet->Det()*6.0);
00248       }
00249     }    
00250     return sum;
00251   }
00252 
00253   // Operator applied to a vector U_array
00254   static inline double Local_matrix_bo_cell(const BoCeData* bocedata,
00255                                           double*  U_array, int num_v) {
00256     static double sum;
00257     Tetraeder_storage* tet;
00258 
00259     sum=0.0;
00260     for(tet=bocedata->Give_tets();tet!=NULL;tet=tet->Next()) {
00261       if(tet->N0()==num_v) {
00262         sum = sum +
00263           XM0*(XM0*UAA0 + XM1*UAA1 + XM2*UAA2 + XM3*UAA3) / (tet->Det()*6.0);
00264       }
00265       else if(tet->N1()==num_v) {
00266         sum = sum +
00267           XM1*(XM0*UAA0 + XM1*UAA1 + XM2*UAA2 + XM3*UAA3) / (tet->Det()*6.0);
00268       }
00269       else if(tet->N2()==num_v) {
00270         sum = sum +
00271           XM2*(XM0*UAA0 + XM1*UAA1 + XM2*UAA2 + XM3*UAA3) / (tet->Det()*6.0);
00272       }
00273       else if(tet->N3()==num_v) {
00274         sum = sum +
00275           XM3*(XM0*UAA0 + XM1*UAA1 + XM2*UAA2 + XM3*UAA3) / (tet->Det()*6.0);
00276       }
00277     }    
00278     return sum;
00279   }
00280 
00281   static inline double diag_F_bo_cell(const BoCeData* bocedata,
00282                                       int num_v) {
00283    static double sum;
00284     Tetraeder_storage* tet;
00285 
00286     sum=0.0;
00287     for(tet=bocedata->Give_tets();tet!=NULL;tet=tet->Next()) {
00288       if(tet->N0()==num_v) {
00289         sum = sum + XM0*XM0 / (tet->Det()*6.0);
00290       }
00291       else if(tet->N1()==num_v) {
00292         sum = sum + XM1*XM1 / (tet->Det()*6.0);
00293       }
00294       else if(tet->N2()==num_v) {
00295         sum = sum + XM2*XM2 / (tet->Det()*6.0);
00296       }
00297       else if(tet->N3()==num_v) {
00298         sum = sum + XM3*XM3 / (tet->Det()*6.0);
00299       }
00300     }     
00301     return sum;
00302   }
00303 };
00304 
00305 
00306 // Formeln fuer dydy
00307 class dydy_FE_const {
00308  public:
00309   diff_type typ_u() { return  Dy_type; }; 
00310   diff_type typ_v() { return  Dy_type; };
00311 
00312   dydy_FE_const() { }
00313   static inline int    ops_interior() { return 18; }
00314   static inline int    ops_diag_interior() { return 1; }
00315   static inline double stencil(double uN,
00316                                double uW, double uM, double uO,
00317                                double uS, double uT, double uD,
00318                                double uND, double uWN, double uWT,
00319                                double uED, double uST, double uES,
00320                                double uEST, double uWND,
00321                                double h_mesh) {
00322     return  (  8.0 * uM - 2.0 * (uN + uS)
00323                - uO  - uW  - uT - uD
00324                - uND - uST - uWN - uES + uWT + uED
00325                + uEST + uWND ) * h_mesh / 3.0; 
00326   }
00327   static inline double diag_stencil(double h_mesh) {
00328     return 8.0 / 3.0 * h_mesh;
00329   }
00330   static inline double F_reg_cell(double*const * u_Recell, int num_var,
00331                                   dir_sons dir_v,
00332                                   double h_mesh, double* sten) {
00333     static int i;
00334     static double sum;
00335 
00336     sum = 0.0;
00337     for(i=0;i<8;++i) sum = sum + sten[dir_v+8*i] * u_Recell[i][num_var];
00338     return sum * h_mesh / 6.0;
00339   }
00340   static inline double Factor_reg_cell(double h_mesh) {
00341     return h_mesh / 6.0;
00342   }
00343   static inline double diag_F_reg_cell(dir_sons dir_v,
00344                                        double h_mesh, double* sten) {
00345     return sten[dir_v*9] * h_mesh / 6.0;
00346   }
00347 
00348   // Operator applied to variable with number
00349   static inline double F_bo_cell(const BoCeData* bocedata,
00350                                  int num_v, int num_variable) {
00351     static double sum;
00352     Tetraeder_storage* tet;
00353 
00354     sum=0.0;
00355     for(tet=bocedata->Give_tets();tet!=NULL;tet=tet->Next()) {
00356       if(tet->N0()==num_v) {
00357         sum = sum +
00358               YM0*(YM0*UU0 + YM1*UU1 + YM2*UU2 + YM3*UU3) / (tet->Det()*6.0);
00359       }
00360       else if(tet->N1()==num_v) {
00361         sum = sum +
00362               YM1*(YM0*UU0 + YM1*UU1 + YM2*UU2 + YM3*UU3) / (tet->Det()*6.0);
00363       }
00364       else if(tet->N2()==num_v) {
00365         sum = sum +
00366               YM2*(YM0*UU0 + YM1*UU1 + YM2*UU2 + YM3*UU3) / (tet->Det()*6.0);
00367       }
00368       else if(tet->N3()==num_v) {
00369         sum = sum +
00370               YM3*(YM0*UU0 + YM1*UU1 + YM2*UU2 + YM3*UU3) / (tet->Det()*6.0);
00371       }
00372     }    
00373     return sum;
00374   }
00375   // Operator applied to a vector U_array
00376   static inline double Local_matrix_bo_cell(const BoCeData* bocedata,
00377                                           double*  U_array, int num_v) {
00378     static double sum;
00379     Tetraeder_storage* tet;
00380 
00381     sum=0.0;
00382     for(tet=bocedata->Give_tets();tet!=NULL;tet=tet->Next()) {
00383       if(tet->N0()==num_v) {
00384         sum = sum +
00385           YM0*(YM0*UAA0 + YM1*UAA1 + YM2*UAA2 + YM3*UAA3) / (tet->Det()*6.0);
00386       }
00387       else if(tet->N1()==num_v) {
00388         sum = sum +
00389           YM1*(YM0*UAA0 + YM1*UAA1 + YM2*UAA2 + YM3*UAA3) / (tet->Det()*6.0);
00390       }
00391       else if(tet->N2()==num_v) {
00392         sum = sum +
00393           YM2*(YM0*UAA0 + YM1*UAA1 + YM2*UAA2 + YM3*UAA3) / (tet->Det()*6.0);
00394       }
00395       else if(tet->N3()==num_v) {
00396         sum = sum +
00397           YM3*(YM0*UAA0 + YM1*UAA1 + YM2*UAA2 + YM3*UAA3) / (tet->Det()*6.0);
00398       }
00399     }    
00400     return sum;
00401   }
00402   static inline double diag_F_bo_cell(const BoCeData* bocedata,
00403                                       int num_v) {
00404    static double sum;
00405     Tetraeder_storage* tet;
00406 
00407     sum=0.0;
00408     for(tet=bocedata->Give_tets();tet!=NULL;tet=tet->Next()) {
00409       if(tet->N0()==num_v) {
00410         sum = sum + YM0*YM0 / (tet->Det()*6.0);
00411       }
00412       else if(tet->N1()==num_v) {
00413         sum = sum + YM1*YM1 / (tet->Det()*6.0);
00414       }
00415       else if(tet->N2()==num_v) {
00416         sum = sum + YM2*YM2 / (tet->Det()*6.0);
00417       }
00418       else if(tet->N3()==num_v) {
00419         sum = sum + YM3*YM3 / (tet->Det()*6.0);
00420       }
00421     }     
00422     return sum;
00423   }
00424 };
00425 
00426 
00427 
00428 // Formeln fuer dzdz
00429 class dzdz_FE_const {
00430  public:  
00431   diff_type typ_u() { return  Dz_type; }; 
00432   diff_type typ_v() { return  Dz_type; };
00433 
00434   dzdz_FE_const() { }
00435   static inline int    ops_interior() { return 4; }
00436   static inline int    ops_diag_interior() { return 1; }
00437   static inline double stencil(double uN,
00438                                double uW, double uM, double uO,
00439                                double uS, double uT, double uD,
00440                                double uND, double uWN, double uWT,
00441                                double uED, double uST, double uES,
00442                                double uEST, double uWND,
00443                                double h_mesh) {
00444     return  (  2.0 * uM - uT - uD) * h_mesh; 
00445   }
00446   static inline double diag_stencil(double h_mesh) {
00447     return 2.0 * h_mesh; 
00448   }
00449   static inline double F_reg_cell(double*const * u_Recell, int num_var,
00450                                   dir_sons dir_v,
00451                                   double h_mesh, double* sten) {
00452     static int i;
00453     static double sum;
00454 
00455     sum = 0.0;
00456     for(i=0;i<8;++i) sum = sum + sten[dir_v+8*i] * u_Recell[i][num_var];
00457     return sum * h_mesh / 6.0;
00458   }
00459   static inline double Factor_reg_cell(double h_mesh) {
00460     return h_mesh / 6.0;
00461   }
00462   static inline double diag_F_reg_cell(dir_sons dir_v,
00463                                        double h_mesh, double* sten) {
00464     return sten[dir_v*9] * h_mesh / 6.0;
00465   }
00466 
00467   // Operator applied to variable with number
00468   static inline double F_bo_cell(const BoCeData* bocedata,
00469                                  int num_v, int num_variable) {
00470     static double sum;
00471     Tetraeder_storage* tet;
00472 
00473     sum=0.0;
00474     for(tet=bocedata->Give_tets();tet!=NULL;tet=tet->Next()) {
00475       if(tet->N0()==num_v) {
00476         sum = sum +
00477               ZM0*(ZM0*UU0 + ZM1*UU1 + ZM2*UU2 + ZM3*UU3) / (tet->Det()*6.0);
00478       }
00479       else if(tet->N1()==num_v) {
00480         sum = sum +
00481               ZM1*(ZM0*UU0 + ZM1*UU1 + ZM2*UU2 + ZM3*UU3) / (tet->Det()*6.0);
00482       }
00483       else if(tet->N2()==num_v) {
00484         sum = sum +
00485               ZM2*(ZM0*UU0 + ZM1*UU1 + ZM2*UU2 + ZM3*UU3) / (tet->Det()*6.0);
00486       }
00487       else if(tet->N3()==num_v) {
00488         sum = sum +
00489               ZM3*(ZM0*UU0 + ZM1*UU1 + ZM2*UU2 + ZM3*UU3) / (tet->Det()*6.0);
00490       }
00491     }     
00492     return sum;
00493   }
00494   // Operator applied to a vector U_array
00495   static inline double Local_matrix_bo_cell(const BoCeData* bocedata,
00496                                           double*  U_array, int num_v) {
00497     static double sum;
00498     Tetraeder_storage* tet;
00499 
00500     sum=0.0;
00501     for(tet=bocedata->Give_tets();tet!=NULL;tet=tet->Next()) {
00502       if(tet->N0()==num_v) {
00503         sum = sum +
00504           ZM0*(ZM0*UAA0 + ZM1*UAA1 + ZM2*UAA2 + ZM3*UAA3) / (tet->Det()*6.0);
00505       }
00506       else if(tet->N1()==num_v) {
00507         sum = sum +
00508           ZM1*(ZM0*UAA0 + ZM1*UAA1 + ZM2*UAA2 + ZM3*UAA3) / (tet->Det()*6.0);
00509       }
00510       else if(tet->N2()==num_v) {
00511         sum = sum +
00512           ZM2*(ZM0*UAA0 + ZM1*UAA1 + ZM2*UAA2 + ZM3*UAA3) / (tet->Det()*6.0);
00513       }
00514       else if(tet->N3()==num_v) {
00515         sum = sum +
00516           ZM3*(ZM0*UAA0 + ZM1*UAA1 + ZM2*UAA2 + ZM3*UAA3) / (tet->Det()*6.0);
00517       }
00518     }     
00519     return sum;
00520   }
00521   static inline double diag_F_bo_cell(const BoCeData* bocedata,
00522                                       int num_v) {
00523    static double sum;
00524     Tetraeder_storage* tet;
00525 
00526     sum=0.0;
00527     for(tet=bocedata->Give_tets();tet!=NULL;tet=tet->Next()) {
00528       if(tet->N0()==num_v) {
00529         sum = sum + ZM0*ZM0 / (tet->Det()*6.0);
00530       }
00531       else if(tet->N1()==num_v) {
00532         sum = sum + ZM1*ZM1 / (tet->Det()*6.0);
00533       }
00534       else if(tet->N2()==num_v) {
00535         sum = sum + ZM2*ZM2 / (tet->Det()*6.0);
00536       }
00537       else if(tet->N3()==num_v) {
00538         sum = sum + ZM3*ZM3 / (tet->Det()*6.0);
00539       }
00540     }     
00541     return sum;
00542   }
00543 };
00544 
00545 
00546 //  Helm
00547 class helm_FE_const {
00548  public:
00549   diff_type typ_u() { return  Helm_type; }; 
00550   diff_type typ_v() { return  Helm_type; };
00551 
00552   helm_FE_const() { }
00553   static inline int    ops_interior() { return 21; }
00554   static inline int    ops_diag_interior() { return 3; }
00555   static inline double stencil(double uN,
00556                                double uW, double uM, double uO,
00557                                double uS, double uT, double uD,
00558                                double uND, double uWN, double uWT,
00559                                double uED, double uST, double uES,
00560                                double uEST, double uWND,
00561                                double h_mesh) {
00562     return  (48.0* uM + 6.0*(uT + uD + uO + uW +
00563                              uES + uWN +
00564                              uND + uST) +
00565              4.0*(uWT  + uED  + uN + uS +
00566                   uEST + uWND)) * h_mesh*h_mesh*h_mesh / 120.0; 
00567   }
00568   static inline double diag_stencil(double h_mesh) {
00569     return h_mesh*h_mesh*h_mesh * 2.0 / 5.0;
00570   }
00571   static inline double F_reg_cell(double*const * u_Recell, int num_var,
00572                                   dir_sons dir_v,
00573                                   double h_mesh, double* sten) {
00574     static int i;
00575     static double sum;
00576   
00577     sum = 0.0;
00578     for(i=0;i<8;++i) sum = sum + sten[dir_v+8*i] * u_Recell[i][num_var];
00579     return sum * h_mesh*h_mesh*h_mesh / 120.0;
00580   }
00581   static inline double Factor_reg_cell(double h_mesh) {
00582     return h_mesh*h_mesh*h_mesh / 120.0;
00583   }
00584   static inline double diag_F_reg_cell(dir_sons dir_v,
00585                                        double h_mesh, double* sten) {
00586     return sten[dir_v*9] * h_mesh*h_mesh*h_mesh / 120.0;;
00587   }
00588 
00589   // Operator applied to variable with number
00590   static inline double F_bo_cell(const BoCeData* bocedata,
00591                                  int num_v, int num_variable) {
00592     static double sum;
00593     Tetraeder_storage* tet;
00594     
00595     sum=0.0;
00596     for(tet=bocedata->Give_tets();tet!=NULL;tet=tet->Next()) {
00597       if(tet->N0()==num_v) {
00598         sum = sum + tet->Det() *
00599               (2.0 * UU0 + UU1 + UU2 + UU3) / 120.0;
00600       }
00601       else if(tet->N1()==num_v) {
00602         sum = sum + tet->Det() *
00603               (2.0 * UU1 + UU0 + UU2 + UU3) / 120.0;
00604       }
00605       else if(tet->N2()==num_v) {
00606         sum = sum + tet->Det() *
00607               (2.0 * UU2 + UU0 + UU1 + UU3) / 120.0;
00608       }
00609       else if(tet->N3()==num_v) {
00610         sum = sum + tet->Det() *
00611               (2.0 * UU3 + UU0 + UU1 + UU2) / 120.0;
00612       }
00613     }  
00614     return sum;
00615   }
00616   // Operator applied to a vector U_array
00617   static inline double Local_matrix_bo_cell(const BoCeData* bocedata,
00618                                           double*  U_array, int num_v) {
00619     static double sum;
00620     Tetraeder_storage* tet;
00621 
00622 
00623     sum=0.0;
00624     for(tet=bocedata->Give_tets();tet!=NULL;tet=tet->Next()) {
00625       if(tet->N0()==num_v) {
00626         sum = sum + tet->Det() *
00627               (2.0 * UAA0 + UAA1 + UAA2 + UAA3) / 120.0;
00628       }
00629       else if(tet->N1()==num_v) {
00630         sum = sum + tet->Det() *
00631               (2.0 * UAA1 + UAA0 + UAA2 + UAA3) / 120.0;
00632       }
00633       else if(tet->N2()==num_v) {
00634         sum = sum + tet->Det() *
00635               (2.0 * UAA2 + UAA0 + UAA1 + UAA3) / 120.0;
00636       }
00637       else if(tet->N3()==num_v) {
00638         sum = sum + tet->Det() *
00639               (2.0 * UAA3 + UAA0 + UAA1 + UAA2) / 120.0;
00640       }
00641     }  
00642     return sum;
00643   }
00644   static inline double diag_F_bo_cell(const BoCeData* bocedata,
00645                                       int num_v) {
00646    static double sum;
00647     Tetraeder_storage* tet;
00648     sum=0.0;
00649     for(tet=bocedata->Give_tets();tet!=NULL;tet=tet->Next()) {
00650       if(tet->N0()==num_v || tet->N1()==num_v ||
00651          tet->N2()==num_v || tet->N3()==num_v)  
00652       sum = sum + tet->Det() / 60.0;
00653     }  
00654     return sum;
00655   }
00656 };
00657 
00659 // Unsymmetrische Terme
00661 
00662 
00663 // Formeln fuer dxdy
00664 class dxdy_FE_const {
00665  public:  
00666   diff_type typ_u() { return  Dx_type; }; 
00667   diff_type typ_v() { return  Dy_type; };
00668 
00669   dxdy_FE_const() { }
00670   static inline int    ops_interior() { return 19; }
00671   static inline int    ops_diag_interior() { return 2; }
00672   static inline double stencil(double uN,
00673                                double uW, double uM, double uO,
00674                                double uS, double uT, double uD,
00675                                double uND, double uWN, double uWT,
00676                                double uED, double uST, double uES,
00677                                double uEST, double uWND,
00678                                double h_mesh) {
00679     return  (  8.0 * uM - 4.0 * (uO + uW) 
00680                         + 2.0 * (uWN + uES - uN - uS)
00681                         + (   uWND + uEST 
00682                             - uND  - uST - uT - uD
00683                             + uWT  + uED)) * h_mesh / 6.0; 
00684   }
00685   static inline double diag_stencil(double h_mesh) {
00686     return 8.0 * h_mesh / 6.0; 
00687   }
00688   static inline double F_reg_cell(double*const * u_Recell, int num_var,
00689                                   dir_sons dir_v,
00690                                   double h_mesh, double* sten) {
00691     static int i;
00692     static double sum;
00693 
00694     sum = 0.0;
00695     for(i=0;i<8;++i) sum = sum + sten[dir_v+8*i] * u_Recell[i][num_var];
00696     return sum * h_mesh / 6.0;
00697   }
00698   static inline double Factor_reg_cell(double h_mesh) {
00699     return h_mesh / 6.0;
00700   }
00701   static inline double diag_F_reg_cell(dir_sons dir_v,
00702                                        double h_mesh, double* sten) {
00703     return sten[dir_v*9] * h_mesh / 6.0;
00704   }
00705 
00706   // Operator applied to variable with number
00707   static inline double F_bo_cell(const BoCeData* bocedata,
00708                                  int num_v, int num_variable) {
00709     static double sum;
00710     Tetraeder_storage* tet;
00711 
00712     sum=0.0;
00713     for(tet=bocedata->Give_tets();tet!=NULL;tet=tet->Next()) {
00714       if(tet->N0()==num_v) {
00715         sum = sum +
00716               YM0*(XM0*UU0 + XM1*UU1 + XM2*UU2 + XM3*UU3) / (tet->Det()*6.0);
00717       }
00718       else if(tet->N1()==num_v) {
00719         sum = sum +
00720               YM1*(XM0*UU0 + XM1*UU1 + XM2*UU2 + XM3*UU3) / (tet->Det()*6.0);
00721       }
00722       else if(tet->N2()==num_v) {
00723         sum = sum +
00724               YM2*(XM0*UU0 + XM1*UU1 + XM2*UU2 + XM3*UU3) / (tet->Det()*6.0);
00725       }
00726       else if(tet->N3()==num_v) {
00727         sum = sum +
00728               YM3*(XM0*UU0 + XM1*UU1 + XM2*UU2 + XM3*UU3) / (tet->Det()*6.0);
00729       }
00730     }     
00731     return sum;
00732   }
00733   // Operator applied to a vector U_array
00734   static inline double Local_matrix_bo_cell(const BoCeData* bocedata,
00735                                           double*  U_array, int num_v) {
00736     static double sum;
00737     Tetraeder_storage* tet;
00738 
00739     sum=0.0;
00740     for(tet=bocedata->Give_tets();tet!=NULL;tet=tet->Next()) {
00741       if(tet->N0()==num_v) {
00742         sum = sum +
00743           YM0*(XM0*UAA0 + XM1*UAA1 + XM2*UAA2 + XM3*UAA3) / (tet->Det()*6.0);
00744       }
00745       else if(tet->N1()==num_v) {
00746         sum = sum +
00747           YM1*(XM0*UAA0 + XM1*UAA1 + XM2*UAA2 + XM3*UAA3) / (tet->Det()*6.0);
00748       }
00749       else if(tet->N2()==num_v) {
00750         sum = sum +
00751           YM2*(XM0*UAA0 + XM1*UAA1 + XM2*UAA2 + XM3*UAA3) / (tet->Det()*6.0);
00752       }
00753       else if(tet->N3()==num_v) {
00754         sum = sum +
00755           YM3*(XM0*UAA0 + XM1*UAA1 + XM2*UAA2 + XM3*UAA3) / (tet->Det()*6.0);
00756       }
00757     }     
00758     return sum;
00759   }
00760   static inline double diag_F_bo_cell(const BoCeData* bocedata,
00761                                       int num_v) {
00762    static double sum;
00763     Tetraeder_storage* tet;
00764 
00765     sum=0.0;
00766     for(tet=bocedata->Give_tets();tet!=NULL;tet=tet->Next()) {
00767       if(tet->N0()==num_v) {
00768         sum = sum + YM0*XM0 / (tet->Det()*6.0);
00769       }
00770       else if(tet->N1()==num_v) {
00771         sum = sum + YM1*XM1 / (tet->Det()*6.0);
00772       }
00773       else if(tet->N2()==num_v) {
00774         sum = sum + YM2*XM2 / (tet->Det()*6.0);
00775       }
00776       else if(tet->N3()==num_v) {
00777         sum = sum + YM3*XM3 / (tet->Det()*6.0);
00778       }
00779     }     
00780     return sum;
00781   }
00782 };
00783 
00784 // Formeln fuer dydx
00785 class dydx_FE_const {
00786  public:  
00787   diff_type typ_u() { return  Dy_type; }; 
00788   diff_type typ_v() { return  Dx_type; };
00789 
00790   dydx_FE_const() { }
00791   static inline int    ops_interior() { return 19; }
00792   static inline int    ops_diag_interior() { return 2; }
00793   static inline double stencil(double uN,
00794                                double uW, double uM, double uO,
00795                                double uS, double uT, double uD,
00796                                double uND, double uWN, double uWT,
00797                                double uED, double uST, double uES,
00798                                double uEST, double uWND,
00799                                double h_mesh) {
00800     return  (  8.0 * uM - 4.0 * (uO + uW) 
00801                         + 2.0 * (uWN + uES - uN - uS)
00802                         + (   uWND + uEST 
00803                             - uND  - uST - uT - uD
00804                             + uWT  + uED)) * h_mesh / 6.0; 
00805   }
00806   static inline double diag_stencil(double h_mesh) {
00807     return 8.0 * h_mesh / 6.0; 
00808   }
00809   static inline double F_reg_cell(double*const * u_Recell, int num_var,
00810                                   dir_sons dir_v,
00811                                   double h_mesh, double* sten) {
00812     static int i;
00813     static double sum;
00814 
00815     sum = 0.0;
00816     for(i=0;i<8;++i) sum = sum + sten[dir_v+8*i] * u_Recell[i][num_var];
00817     return sum * h_mesh / 6.0;
00818   }
00819   static inline double Factor_reg_cell(double h_mesh) {
00820     return h_mesh / 6.0;
00821   }
00822   static inline double diag_F_reg_cell(dir_sons dir_v,
00823                                        double h_mesh, double* sten) {
00824     return sten[dir_v*9] * h_mesh / 6.0;
00825   }
00826 
00827   // Operator applied to variable with number
00828   static inline double F_bo_cell(const BoCeData* bocedata,
00829                                  int num_v, int num_variable) {
00830     static double sum;
00831     Tetraeder_storage* tet;
00832 
00833     sum=0.0;
00834     for(tet=bocedata->Give_tets();tet!=NULL;tet=tet->Next()) {
00835       if(tet->N0()==num_v) {
00836         sum = sum +
00837               XM0*(YM0*UU0 + YM1*UU1 + YM2*UU2 + YM3*UU3) / (tet->Det()*6.0);
00838       }
00839       else if(tet->N1()==num_v) {
00840         sum = sum +
00841               XM1*(YM0*UU0 + YM1*UU1 + YM2*UU2 + YM3*UU3) / (tet->Det()*6.0);
00842       }
00843       else if(tet->N2()==num_v) {
00844         sum = sum +
00845               XM2*(YM0*UU0 + YM1*UU1 + YM2*UU2 + YM3*UU3) / (tet->Det()*6.0);
00846       }
00847       else if(tet->N3()==num_v) {
00848         sum = sum +
00849               XM3*(YM0*UU0 + YM1*UU1 + YM2*UU2 + YM3*UU3) / (tet->Det()*6.0);
00850       }
00851     }     
00852     return sum;
00853   }
00854   // Operator applied to a vector U_array
00855   static inline double Local_matrix_bo_cell(const BoCeData* bocedata,
00856                                           double*  U_array, int num_v) {
00857     static double sum;
00858     Tetraeder_storage* tet;
00859 
00860     sum=0.0;
00861     for(tet=bocedata->Give_tets();tet!=NULL;tet=tet->Next()) {
00862       if(tet->N0()==num_v) {
00863         sum = sum +
00864           XM0*(YM0*UAA0 + YM1*UAA1 + YM2*UAA2 + YM3*UAA3) / (tet->Det()*6.0);
00865       }
00866       else if(tet->N1()==num_v) {
00867         sum = sum +
00868           XM1*(YM0*UAA0 + YM1*UAA1 + YM2*UAA2 + YM3*UAA3) / (tet->Det()*6.0);
00869       }
00870       else if(tet->N2()==num_v) {
00871         sum = sum +
00872           XM2*(YM0*UAA0 + YM1*UAA1 + YM2*UAA2 + YM3*UAA3) / (tet->Det()*6.0);
00873       }
00874       else if(tet->N3()==num_v) {
00875         sum = sum +
00876           XM3*(YM0*UAA0 + YM1*UAA1 + YM2*UAA2 + YM3*UAA3) / (tet->Det()*6.0);
00877       }
00878     }     
00879     return sum;
00880   }
00881   static inline double diag_F_bo_cell(const BoCeData* bocedata,
00882                                       int num_v) {
00883    static double sum;
00884     Tetraeder_storage* tet;
00885 
00886     sum=0.0;
00887     for(tet=bocedata->Give_tets();tet!=NULL;tet=tet->Next()) {
00888       if(tet->N0()==num_v) {
00889         sum = sum + XM0*YM0 / (tet->Det()*6.0);
00890       }
00891       else if(tet->N1()==num_v) {
00892         sum = sum + XM1*YM1 / (tet->Det()*6.0);
00893       }
00894       else if(tet->N2()==num_v) {
00895         sum = sum + XM2*YM2 / (tet->Det()*6.0);
00896       }
00897       else if(tet->N3()==num_v) {
00898         sum = sum + XM3*YM3 / (tet->Det()*6.0);
00899       }
00900     }     
00901     return sum;
00902   }
00903 };
00904 
00906 
00907 
00908 
00909 // Formeln fuer dxdz
00910 class dxdz_FE_const {
00911  public:  
00912   diff_type typ_u() { return  Dx_type; }; 
00913   diff_type typ_v() { return  Dz_type; };
00914 
00915   dxdz_FE_const() { }
00916   static inline int    ops_interior() { return 18; }
00917   static inline int    ops_diag_interior() { return 2; }
00918   static inline double stencil(double uN,
00919                                double uW, double uM, double uO,
00920                                double uS, double uT, double uD,
00921                                double uND, double uWN, double uWT,
00922                                double uED, double uST, double uES,
00923                                double uEST, double uWND,
00924                                double h_mesh) {
00925     return  (  4.0 * uM + 2.0 * (uWT + uED - uW - uO - uT - uD)
00926                         + (   uND  + uST + uWN  + uES 
00927                             - uWND - uEST - uN - uS)) * h_mesh / 6.0; 
00928   }
00929   static inline double diag_stencil(double h_mesh) {
00930     return 4.0 * h_mesh / 6.0; 
00931   }
00932   static inline double F_reg_cell(double*const * u_Recell, int num_var,
00933                                   dir_sons dir_v,
00934                                   double h_mesh, double* sten) {
00935     static int i;
00936     static double sum;
00937 
00938     sum = 0.0;
00939     for(i=0;i<8;++i) sum = sum + sten[dir_v+8*i] * u_Recell[i][num_var];
00940     return sum * h_mesh / 6.0;
00941   }
00942   static inline double Factor_reg_cell(double h_mesh) {
00943     return h_mesh / 6.0;
00944   }
00945   static inline double diag_F_reg_cell(dir_sons dir_v,
00946                                        double h_mesh, double* sten) {
00947     return sten[dir_v*9] * h_mesh / 6.0;
00948   }
00949 
00950   // Operator applied to variable with number
00951   static inline double F_bo_cell(const BoCeData* bocedata,
00952                                  int num_v, int num_variable) {
00953     static double sum;
00954     Tetraeder_storage* tet;
00955 
00956     sum=0.0;
00957     for(tet=bocedata->Give_tets();tet!=NULL;tet=tet->Next()) {
00958       if(tet->N0()==num_v) {
00959         sum = sum +
00960               ZM0*(XM0*UU0 + XM1*UU1 + XM2*UU2 + XM3*UU3) / (tet->Det()*6.0);
00961       }
00962       else if(tet->N1()==num_v) {
00963         sum = sum +
00964               ZM1*(XM0*UU0 + XM1*UU1 + XM2*UU2 + XM3*UU3) / (tet->Det()*6.0);
00965       }
00966       else if(tet->N2()==num_v) {
00967         sum = sum +
00968               ZM2*(XM0*UU0 + XM1*UU1 + XM2*UU2 + XM3*UU3) / (tet->Det()*6.0);
00969       }
00970       else if(tet->N3()==num_v) {
00971         sum = sum +
00972               ZM3*(XM0*UU0 + XM1*UU1 + XM2*UU2 + XM3*UU3) / (tet->Det()*6.0);
00973       }
00974     }     
00975     return sum;
00976   }
00977   // Operator applied to a vector U_array
00978   static inline double Local_matrix_bo_cell(const BoCeData* bocedata,
00979                                           double*  U_array, int num_v) {
00980     static double sum;
00981     Tetraeder_storage* tet;
00982 
00983     sum=0.0;
00984     for(tet=bocedata->Give_tets();tet!=NULL;tet=tet->Next()) {
00985       if(tet->N0()==num_v) {
00986         sum = sum +
00987           ZM0*(XM0*UAA0 + XM1*UAA1 + XM2*UAA2 + XM3*UAA3) / (tet->Det()*6.0);
00988       }
00989       else if(tet->N1()==num_v) {
00990         sum = sum +
00991           ZM1*(XM0*UAA0 + XM1*UAA1 + XM2*UAA2 + XM3*UAA3) / (tet->Det()*6.0);
00992       }
00993       else if(tet->N2()==num_v) {
00994         sum = sum +
00995           ZM2*(XM0*UAA0 + XM1*UAA1 + XM2*UAA2 + XM3*UAA3) / (tet->Det()*6.0);
00996       }
00997       else if(tet->N3()==num_v) {
00998         sum = sum +
00999           ZM3*(XM0*UAA0 + XM1*UAA1 + XM2*UAA2 + XM3*UAA3) / (tet->Det()*6.0);
01000       }
01001     }     
01002     return sum;
01003   }
01004   static inline double diag_F_bo_cell(const BoCeData* bocedata,
01005                                       int num_v) {
01006    static double sum;
01007     Tetraeder_storage* tet;
01008 
01009     sum=0.0;
01010     for(tet=bocedata->Give_tets();tet!=NULL;tet=tet->Next()) {
01011       if(tet->N0()==num_v) {
01012         sum = sum + ZM0*XM0 / (tet->Det()*6.0);
01013       }
01014       else if(tet->N1()==num_v) {
01015         sum = sum + ZM1*XM1 / (tet->Det()*6.0);
01016       }
01017       else if(tet->N2()==num_v) {
01018         sum = sum + ZM2*XM2 / (tet->Det()*6.0);
01019       }
01020       else if(tet->N3()==num_v) {
01021         sum = sum + ZM3*XM3 / (tet->Det()*6.0);
01022       }
01023     }     
01024     return sum;
01025   }
01026 };
01027 
01028 // Formeln fuer dzdx
01029 class dzdx_FE_const {
01030  public:  
01031   diff_type typ_u() { return  Dz_type; }; 
01032   diff_type typ_v() { return  Dx_type; };
01033 
01034   dzdx_FE_const() { }
01035   static inline int    ops_interior() { return 18; }
01036   static inline int    ops_diag_interior() { return 2; }
01037   static inline double stencil(double uN,
01038                                double uW, double uM, double uO,
01039                                double uS, double uT, double uD,
01040                                double uND, double uWN, double uWT,
01041                                double uED, double uST, double uES,
01042                                double uEST, double uWND,
01043                                double h_mesh) {
01044     return  (  4.0 * uM + 2.0 * (uWT + uED - uW - uO - uT - uD)
01045                         + (   uND  + uST + uWN  + uES 
01046                             - uWND - uEST - uN - uS)) * h_mesh / 6.0; 
01047   }
01048   static inline double diag_stencil(double h_mesh) {
01049     return 4.0 * h_mesh / 6.0; 
01050   }
01051   static inline double F_reg_cell(double*const * u_Recell, int num_var,
01052                                   dir_sons dir_v,
01053                                   double h_mesh, double* sten) {
01054     static int i;
01055     static double sum;
01056 
01057     sum = 0.0;
01058     for(i=0;i<8;++i) sum = sum + sten[dir_v+8*i] * u_Recell[i][num_var];
01059     return sum * h_mesh / 6.0;
01060   }
01061   static inline double Factor_reg_cell(double h_mesh) {
01062     return h_mesh / 6.0;
01063   }
01064   static inline double diag_F_reg_cell(dir_sons dir_v,
01065                                        double h_mesh, double* sten) {
01066     return sten[dir_v*9] * h_mesh / 6.0;
01067   }
01068 
01069   // Operator applied to variable with number
01070   static inline double F_bo_cell(const BoCeData* bocedata,
01071                                  int num_v, int num_variable) {
01072     static double sum;
01073     Tetraeder_storage* tet;
01074 
01075     sum=0.0;
01076     for(tet=bocedata->Give_tets();tet!=NULL;tet=tet->Next()) {
01077       if(tet->N0()==num_v) {
01078         sum = sum +
01079               XM0*(ZM0*UU0 + ZM1*UU1 + ZM2*UU2 + ZM3*UU3) / (tet->Det()*6.0);
01080       }
01081       else if(tet->N1()==num_v) {
01082         sum = sum +
01083               XM1*(ZM0*UU0 + ZM1*UU1 + ZM2*UU2 + ZM3*UU3) / (tet->Det()*6.0);
01084       }
01085       else if(tet->N2()==num_v) {
01086         sum = sum +
01087               XM2*(ZM0*UU0 + ZM1*UU1 + ZM2*UU2 + ZM3*UU3) / (tet->Det()*6.0);
01088       }
01089       else if(tet->N3()==num_v) {
01090         sum = sum +
01091               XM3*(ZM0*UU0 + ZM1*UU1 + ZM2*UU2 + ZM3*UU3) / (tet->Det()*6.0);
01092       }
01093     }     
01094     return sum;
01095   }
01096   // Operator applied to a vector U_array
01097   static inline double Local_matrix_bo_cell(const BoCeData* bocedata,
01098                                           double*  U_array, int num_v) {
01099     static double sum;
01100     Tetraeder_storage* tet;
01101 
01102     sum=0.0;
01103     for(tet=bocedata->Give_tets();tet!=NULL;tet=tet->Next()) {
01104       if(tet->N0()==num_v) {
01105         sum = sum +
01106           XM0*(ZM0*UAA0 + ZM1*UAA1 + ZM2*UAA2 + ZM3*UAA3) / (tet->Det()*6.0);
01107       }
01108       else if(tet->N1()==num_v) {
01109         sum = sum +
01110           XM1*(ZM0*UAA0 + ZM1*UAA1 + ZM2*UAA2 + ZM3*UAA3) / (tet->Det()*6.0);
01111       }
01112       else if(tet->N2()==num_v) {
01113         sum = sum +
01114           XM2*(ZM0*UAA0 + ZM1*UAA1 + ZM2*UAA2 + ZM3*UAA3) / (tet->Det()*6.0);
01115       }
01116       else if(tet->N3()==num_v) {
01117         sum = sum +
01118           XM3*(ZM0*UAA0 + ZM1*UAA1 + ZM2*UAA2 + ZM3*UAA3) / (tet->Det()*6.0);
01119       }
01120     }     
01121     return sum;
01122   }
01123   static inline double diag_F_bo_cell(const BoCeData* bocedata,
01124                                       int num_v) {
01125    static double sum;
01126     Tetraeder_storage* tet;
01127 
01128     sum=0.0;
01129     for(tet=bocedata->Give_tets();tet!=NULL;tet=tet->Next()) {
01130       if(tet->N0()==num_v) {
01131         sum = sum + XM0*ZM0 / (tet->Det()*6.0);
01132       }
01133       else if(tet->N1()==num_v) {
01134         sum = sum + XM1*ZM1 / (tet->Det()*6.0);
01135       }
01136       else if(tet->N2()==num_v) {
01137         sum = sum + XM2*ZM2 / (tet->Det()*6.0);
01138       }
01139       else if(tet->N3()==num_v) {
01140         sum = sum + XM3*ZM3 / (tet->Det()*6.0);
01141       }
01142     }     
01143     return sum;
01144   }
01145 };
01146 
01147 
01149 
01150 
01151 // Formeln fuer dydz
01152 class dydz_FE_const {
01153  public:  
01154   diff_type typ_u() { return  Dy_type; }; 
01155   diff_type typ_v() { return  Dz_type; };
01156 
01157   dydz_FE_const() { }
01158   static inline int    ops_interior() { return 19; }
01159   static inline int    ops_diag_interior() { return 2; }
01160   static inline double stencil(double uN,
01161                                double uW, double uM, double uO,
01162                                double uS, double uT, double uD,
01163                                double uND, double uWN, double uWT,
01164                                double uED, double uST, double uES,
01165                                double uEST, double uWND,
01166                                double h_mesh) {
01167     return  (  8.0 * uM - 4.0 * (uT  + uD)
01168                         + 2.0 * (uND + uST - uN - uS)
01169                         + (   uWND + uEST + uWT + uED
01170                              - uWN - uES  - uW  - uO  )) * h_mesh / 6.0; 
01171   }
01172   static inline double diag_stencil(double h_mesh) {
01173     return 8.0 * h_mesh / 6.0; 
01174   }
01175   static inline double F_reg_cell(double*const * u_Recell, int num_var,
01176                                   dir_sons dir_v,
01177                                   double h_mesh, double* sten) {
01178     static int i;
01179     static double sum;
01180 
01181     sum = 0.0;
01182     for(i=0;i<8;++i) sum = sum + sten[dir_v+8*i] * u_Recell[i][num_var];
01183     return sum * h_mesh / 6.0;
01184   }
01185   static inline double Factor_reg_cell(double h_mesh) {
01186     return h_mesh / 6.0;
01187   }
01188   static inline double diag_F_reg_cell(dir_sons dir_v,
01189                                        double h_mesh, double* sten) {
01190     return sten[dir_v*9] * h_mesh / 6.0;
01191   }
01192 
01193   // Operator applied to variable with number
01194   static inline double F_bo_cell(const BoCeData* bocedata,
01195                                  int num_v, int num_variable) {
01196     static double sum;
01197     Tetraeder_storage* tet;
01198 
01199     sum=0.0;
01200     for(tet=bocedata->Give_tets();tet!=NULL;tet=tet->Next()) {
01201       if(tet->N0()==num_v) {
01202         sum = sum +
01203               ZM0*(YM0*UU0 + YM1*UU1 + YM2*UU2 + YM3*UU3) / (tet->Det()*6.0);
01204       }
01205       else if(tet->N1()==num_v) {
01206         sum = sum +
01207               ZM1*(YM0*UU0 + YM1*UU1 + YM2*UU2 + YM3*UU3) / (tet->Det()*6.0);
01208       }
01209       else if(tet->N2()==num_v) {
01210         sum = sum +
01211               ZM2*(YM0*UU0 + YM1*UU1 + YM2*UU2 + YM3*UU3) / (tet->Det()*6.0);
01212       }
01213       else if(tet->N3()==num_v) {
01214         sum = sum +
01215               ZM3*(YM0*UU0 + YM1*UU1 + YM2*UU2 + YM3*UU3) / (tet->Det()*6.0);
01216       }
01217     }     
01218     return sum;
01219   }
01220   // Operator applied to a vector U_array
01221   static inline double Local_matrix_bo_cell(const BoCeData* bocedata,
01222                                           double*  U_array, int num_v) {
01223     static double sum;
01224     Tetraeder_storage* tet;
01225 
01226     sum=0.0;
01227     for(tet=bocedata->Give_tets();tet!=NULL;tet=tet->Next()) {
01228       if(tet->N0()==num_v) {
01229         sum = sum +
01230           ZM0*(YM0*UAA0 + YM1*UAA1 + YM2*UAA2 + YM3*UAA3) / (tet->Det()*6.0);
01231       }
01232       else if(tet->N1()==num_v) {
01233         sum = sum +
01234           ZM1*(YM0*UAA0 + YM1*UAA1 + YM2*UAA2 + YM3*UAA3) / (tet->Det()*6.0);
01235       }
01236       else if(tet->N2()==num_v) {
01237         sum = sum +
01238           ZM2*(YM0*UAA0 + YM1*UAA1 + YM2*UAA2 + YM3*UAA3) / (tet->Det()*6.0);
01239       }
01240       else if(tet->N3()==num_v) {
01241         sum = sum +
01242           ZM3*(YM0*UAA0 + YM1*UAA1 + YM2*UAA2 + YM3*UAA3) / (tet->Det()*6.0);
01243       }
01244     }     
01245     return sum;
01246   }
01247   static inline double diag_F_bo_cell(const BoCeData* bocedata,
01248                                       int num_v) {
01249    static double sum;
01250     Tetraeder_storage* tet;
01251 
01252     sum=0.0;
01253     for(tet=bocedata->Give_tets();tet!=NULL;tet=tet->Next()) {
01254       if(tet->N0()==num_v) {
01255         sum = sum + ZM0*YM0 / (tet->Det()*6.0);
01256       }
01257       else if(tet->N1()==num_v) {
01258         sum = sum + ZM1*YM1 / (tet->Det()*6.0);
01259       }
01260       else if(tet->N2()==num_v) {
01261         sum = sum + ZM2*YM2 / (tet->Det()*6.0);
01262       }
01263       else if(tet->N3()==num_v) {
01264         sum = sum + ZM3*YM3 / (tet->Det()*6.0);
01265       }
01266     }     
01267     return sum;
01268   }
01269 };
01270 
01271 // Formeln fuer dzdy
01272 class dzdy_FE_const {
01273  public:  
01274   diff_type typ_u() { return  Dz_type; }; 
01275   diff_type typ_v() { return  Dy_type; };
01276 
01277   dzdy_FE_const() { }
01278   static inline int    ops_interior() { return 19; }
01279   static inline int    ops_diag_interior() { return 2; }
01280   static inline double stencil(double uN,
01281                                double uW, double uM, double uO,
01282                                double uS, double uT, double uD,
01283                                double uND, double uWN, double uWT,
01284                                double uED, double uST, double uES,
01285                                double uEST, double uWND,
01286                                double h_mesh) {
01287    return  (  8.0 * uM - 4.0 * (uT  + uD)
01288                         + 2.0 * (uND + uST - uN - uS)
01289                         + (   uWND + uEST + uWT + uED
01290                              - uWN - uES  - uW  - uO  )) * h_mesh / 6.0; 
01291   }
01292   static inline double diag_stencil(double h_mesh) {
01293     return 8.0 * h_mesh / 6.0; 
01294   }
01295   static inline double F_reg_cell(double*const * u_Recell, int num_var,
01296                                   dir_sons dir_v,
01297                                   double h_mesh, double* sten) {
01298     static int i;
01299     static double sum;
01300 
01301     sum = 0.0;
01302     for(i=0;i<8;++i) sum = sum + sten[dir_v+8*i] * u_Recell[i][num_var];
01303     return sum * h_mesh / 6.0;
01304   }
01305   static inline double Factor_reg_cell(double h_mesh) {
01306     return h_mesh / 6.0;
01307   }
01308   static inline double diag_F_reg_cell(dir_sons dir_v,
01309                                        double h_mesh, double* sten) {
01310     return sten[dir_v*9] * h_mesh / 6.0;
01311   }
01312 
01313   // Operator applied to variable with number
01314   static inline double F_bo_cell(const BoCeData* bocedata,
01315                                  int num_v, int num_variable) {
01316     static double sum;
01317     Tetraeder_storage* tet;
01318 
01319     sum=0.0;
01320     for(tet=bocedata->Give_tets();tet!=NULL;tet=tet->Next()) {
01321       if(tet->N0()==num_v) {
01322         sum = sum +
01323               YM0*(ZM0*UU0 + ZM1*UU1 + ZM2*UU2 + ZM3*UU3) / (tet->Det()*6.0);
01324       }
01325       else if(tet->N1()==num_v) {
01326         sum = sum +
01327               YM1*(ZM0*UU0 + ZM1*UU1 + ZM2*UU2 + ZM3*UU3) / (tet->Det()*6.0);
01328       }
01329       else if(tet->N2()==num_v) {
01330         sum = sum +
01331               YM2*(ZM0*UU0 + ZM1*UU1 + ZM2*UU2 + ZM3*UU3) / (tet->Det()*6.0);
01332       }
01333       else if(tet->N3()==num_v) {
01334         sum = sum +
01335               YM3*(ZM0*UU0 + ZM1*UU1 + ZM2*UU2 + ZM3*UU3) / (tet->Det()*6.0);
01336       }
01337     }     
01338     return sum;
01339   }
01340   // Operator applied to a vector U_array
01341   static inline double Local_matrix_bo_cell(const BoCeData* bocedata,
01342                                           double*  U_array, int num_v) {
01343     static double sum;
01344     Tetraeder_storage* tet;
01345 
01346     sum=0.0;
01347     for(tet=bocedata->Give_tets();tet!=NULL;tet=tet->Next()) {
01348       if(tet->N0()==num_v) {
01349         sum = sum +
01350           YM0*(ZM0*UAA0 + ZM1*UAA1 + ZM2*UAA2 + ZM3*UAA3) / (tet->Det()*6.0);
01351       }
01352       else if(tet->N1()==num_v) {
01353         sum = sum +
01354           YM1*(ZM0*UAA0 + ZM1*UAA1 + ZM2*UAA2 + ZM3*UAA3) / (tet->Det()*6.0);
01355       }
01356       else if(tet->N2()==num_v) {
01357         sum = sum +
01358           YM2*(ZM0*UAA0 + ZM1*UAA1 + ZM2*UAA2 + ZM3*UAA3) / (tet->Det()*6.0);
01359       }
01360       else if(tet->N3()==num_v) {
01361         sum = sum +
01362           YM3*(ZM0*UAA0 + ZM1*UAA1 + ZM2*UAA2 + ZM3*UAA3) / (tet->Det()*6.0);
01363       }
01364     }     
01365     return sum;
01366   }
01367   static inline double diag_F_bo_cell(const BoCeData* bocedata,
01368                                       int num_v) {
01369    static double sum;
01370     Tetraeder_storage* tet;
01371 
01372     sum=0.0;
01373     for(tet=bocedata->Give_tets();tet!=NULL;tet=tet->Next()) {
01374       if(tet->N0()==num_v) {
01375         sum = sum + YM0*ZM0 / (tet->Det()*6.0);
01376       }
01377       else if(tet->N1()==num_v) {
01378         sum = sum + YM1*ZM1 / (tet->Det()*6.0);
01379       }
01380       else if(tet->N2()==num_v) {
01381         sum = sum + YM2*ZM2 / (tet->Det()*6.0);
01382       }
01383       else if(tet->N3()==num_v) {
01384         sum = sum + YM3*ZM3 / (tet->Det()*6.0);
01385       }
01386     }     
01387     return sum;
01388   }
01389 };
01390 
01392 
01393 // Formeln fuer dxhelm
01394 class dxhelm_FE_const {
01395  public:
01396   diff_type typ_u() { return  Dx_type; }; 
01397   diff_type typ_v() { return  Helm_type; };
01398 
01399   dxhelm_FE_const() { }
01400   static inline int    ops_interior() { return 17; }
01401   static inline int    ops_diag_interior() { return 0; }
01402   static inline double stencil(double uN,
01403                                double uW, double uM, double uO,
01404                                double uS, double uT, double uD,
01405                                double uND, double uWN, double uWT,
01406                                double uED, double uST, double uES,
01407                                double uEST, double uWND,
01408                                double h_mesh) {
01409     //   return 0.0;
01410     return  (3.0*(uO - uW) 
01411              + uN  - uWN  + uES  - uS
01412              + uND - uWND + uEST - uST
01413              + uED - uD   + uT   - uWT)
01414       * h_mesh * h_mesh / 12.0; 
01415   }
01416   static inline double diag_stencil(double h_mesh) {
01417     return 0.0;
01418   }
01419   static inline double F_reg_cell(double*const * u_Recell, int num_var,
01420                                   dir_sons dir_v,
01421                                   double h_mesh, double* sten) {
01422     static int i;
01423     static double sum;
01424   
01425     sum = 0.0;
01426     for(i=0;i<8;++i) sum = sum + sten[dir_v+8*i] * u_Recell[i][num_var];
01427     // return 0.0;
01428     return sum * h_mesh * h_mesh / 24.0;
01429   }
01430   static inline double Factor_reg_cell(double h_mesh) {
01431     return h_mesh * h_mesh / 24.0;
01432   }
01433   static inline double diag_F_reg_cell(dir_sons dir_v,
01434                                        double h_mesh, double* sten) {
01435     return sten[dir_v*9] * h_mesh * h_mesh / 24.0;
01436   }
01437 
01438   // Operator applied to variable with number
01439   static inline double F_bo_cell(const BoCeData* bocedata,
01440                                  int num_v, int num_variable) {
01441     static double sum;
01442     Tetraeder_storage* tet;
01443 
01444     sum=0.0;
01445     for(tet=bocedata->Give_tets();tet!=NULL;tet=tet->Next()) {
01446       if(tet->N0()==num_v || tet->N1()==num_v  ||
01447          tet->N2()==num_v || tet->N3()==num_v ) {
01448         sum = sum +
01449               (XM0*UU0 + XM1*UU1 + XM2*UU2 + XM3*UU3) / 24.0;
01450       }
01451     }
01452     return sum;
01453   }
01454   // Operator applied to a vector U_array
01455   static inline double Local_matrix_bo_cell(const BoCeData* bocedata,
01456                                           double*  U_array, int num_v) {
01457     static double sum;
01458     Tetraeder_storage* tet;
01459 
01460     sum=0.0;
01461     for(tet=bocedata->Give_tets();tet!=NULL;tet=tet->Next()) {
01462       if(tet->N0()==num_v || tet->N1()==num_v  ||
01463          tet->N2()==num_v || tet->N3()==num_v ) {
01464         sum = sum +
01465               (XM0*UAA0 + XM1*UAA1 + XM2*UAA2 + XM3*UAA3) / 24.0;
01466       }
01467     }
01468     return sum;
01469   }
01470   static inline double diag_F_bo_cell(const BoCeData* bocedata,
01471                                       int num_v) {
01472    static double sum;
01473     Tetraeder_storage* tet;
01474 
01475     sum=0.0;
01476     for(tet=bocedata->Give_tets();tet!=NULL;tet=tet->Next()) {
01477       if(tet->N0()==num_v) sum = sum + XM0 / 24.0;
01478       if(tet->N1()==num_v) sum = sum + XM1 / 24.0;
01479       if(tet->N2()==num_v) sum = sum + XM2 / 24.0;
01480       if(tet->N3()==num_v) sum = sum + XM3 / 24.0;
01481     }
01482     return sum;
01483   }
01484 };
01485 
01486 // Formeln fuer helmdx
01487 class helmdx_FE_const {
01488  public:
01489   diff_type typ_v() { return  Dx_type; }; 
01490   diff_type typ_u() { return  Helm_type; };
01491 
01492   helmdx_FE_const() { }
01493   static inline int    ops_interior() { return 17; }
01494   static inline int    ops_diag_interior() { return 0; }
01495   static inline double stencil(double uN,
01496                                double uW, double uM, double uO,
01497                                double uS, double uT, double uD,
01498                                double uND, double uWN, double uWT,
01499                                double uED, double uST, double uES,
01500                                double uEST, double uWND,
01501                                double h_mesh) {
01502     return  -(3.0*(uO - uW) 
01503               + uN  - uWN  + uES  - uS
01504               + uND - uWND + uEST - uST
01505               + uED - uD   + uT   - uWT)
01506       * h_mesh * h_mesh / 12.0; 
01507   }
01508   static inline double diag_stencil(double h_mesh) {
01509     return 0.0;
01510   }
01511   static inline double F_reg_cell(double*const * u_Recell, int num_var,
01512                                   dir_sons dir_v,
01513                                   double h_mesh, double* sten) {
01514     static int i;
01515     static double sum;
01516   
01517     sum = 0.0;
01518     for(i=0;i<8;++i) sum = sum + sten[dir_v+8*i] * u_Recell[i][num_var];
01519     return sum * h_mesh * h_mesh / 24.0;
01520   }
01521   static inline double Factor_reg_cell(double h_mesh) {
01522     return h_mesh * h_mesh / 24.0;
01523   }
01524   static inline double diag_F_reg_cell(dir_sons dir_v,
01525                                        double h_mesh, double* sten) {
01526     return sten[dir_v*9] * h_mesh * h_mesh / 24.0;
01527   }
01528 
01529   // Operator applied to variable with number
01530   static inline double F_bo_cell(const BoCeData* bocedata,
01531                                  int num_v, int num_variable) {
01532     static double sum;
01533     Tetraeder_storage* tet;
01534 
01535     sum=0.0;
01536     
01537     for(tet=bocedata->Give_tets();tet!=NULL;tet=tet->Next()) {
01538       if(tet->N0()==num_v) sum += XM0*(UU0 + UU1 + UU2 + UU3) / 24.0;
01539       if(tet->N1()==num_v) sum += XM1*(UU0 + UU1 + UU2 + UU3) / 24.0;
01540       if(tet->N2()==num_v) sum += XM2*(UU0 + UU1 + UU2 + UU3) / 24.0;
01541       if(tet->N3()==num_v) sum += XM3*(UU0 + UU1 + UU2 + UU3) / 24.0;
01542     }
01543     return sum;
01544   }
01545   // Operator applied to a vector U_array
01546   static inline double Local_matrix_bo_cell(const BoCeData* bocedata,
01547                                           double*  U_array, int num_v) {
01548    static double sum;
01549     Tetraeder_storage* tet;
01550 
01551     sum=0.0;
01552     
01553     for(tet=bocedata->Give_tets();tet!=NULL;tet=tet->Next()) {
01554       if(tet->N0()==num_v) sum += XM0*(UAA0 + UAA1 + UAA2 + UAA3) / 24.0;
01555       if(tet->N1()==num_v) sum += XM1*(UAA0 + UAA1 + UAA2 + UAA3) / 24.0;
01556       if(tet->N2()==num_v) sum += XM2*(UAA0 + UAA1 + UAA2 + UAA3) / 24.0;
01557       if(tet->N3()==num_v) sum += XM3*(UAA0 + UAA1 + UAA2 + UAA3) / 24.0;
01558     }
01559     return sum;
01560   }
01561   static inline double diag_F_bo_cell(const BoCeData* bocedata,
01562                                       int num_v) {
01563    static double sum;
01564     Tetraeder_storage* tet;
01565 
01566     sum=0.0;
01567     
01568     for(tet=bocedata->Give_tets();tet!=NULL;tet=tet->Next()) {
01569       if(tet->N0()==num_v) sum += XM0 / 24.0;
01570       if(tet->N1()==num_v) sum += XM1 / 24.0;
01571       if(tet->N2()==num_v) sum += XM2 / 24.0;
01572       if(tet->N3()==num_v) sum += XM3 / 24.0;
01573     }
01574     return sum;
01575   }
01576 };
01577 
01578 // Formeln fuer dyhelm
01579 class dyhelm_FE_const {
01580  public:
01581   diff_type typ_u() { return  Dy_type; }; 
01582   diff_type typ_v() { return  Helm_type; };
01583 
01584   dyhelm_FE_const() { }
01585   static inline int    ops_interior() { return 12; }
01586   static inline int    ops_diag_interior() { return 0; }
01587   static inline double stencil(double uN,
01588                                double uW, double uM, double uE,
01589                                double uS, double uT, double uD,
01590                                double uND, double uWN, double uWT,
01591                                double uED, double uST, double uES,
01592                                double uEST, double uWND,
01593                                double h_mesh) {
01594     //    return 0.0;
01595     return  (  uN + uWN + uND + uT + uE 
01596              - uS - uES - uST - uD - uW)
01597       * h_mesh * h_mesh / 6.0; 
01598   }
01599   static inline double diag_stencil(double h_mesh) {
01600     return 0.0;
01601   }
01602   static inline double F_reg_cell(double*const * u_Recell, int num_var,
01603                                   dir_sons dir_v,
01604                                   double h_mesh, double* sten) {
01605     static int i;
01606     static double sum;
01607   
01608     sum = 0.0;
01609     for(i=0;i<8;++i) sum = sum + sten[dir_v+8*i] * u_Recell[i][num_var];
01610     //    return 0.0;
01611     return sum * h_mesh * h_mesh / 24.0;
01612   }
01613   static inline double Factor_reg_cell(double h_mesh) {
01614     return h_mesh * h_mesh / 24.0;
01615   }
01616   static inline double diag_F_reg_cell(dir_sons dir_v,
01617                                        double h_mesh, double* sten) {
01618     return sten[dir_v*9] * h_mesh * h_mesh / 24.0;
01619   }
01620 
01621   // Operator applied to variable with number
01622   static inline double F_bo_cell(const BoCeData* bocedata,
01623                                  int num_v, int num_variable) {
01624     static double sum;
01625     Tetraeder_storage* tet;
01626 
01627     sum=0.0;
01628     for(tet=bocedata->Give_tets();tet!=NULL;tet=tet->Next()) {
01629       if(tet->N0()==num_v || tet->N1()==num_v  ||
01630          tet->N2()==num_v || tet->N3()==num_v ) {
01631         sum = sum +
01632               (YM0*UU0 + YM1*UU1 + YM2*UU2 + YM3*UU3) / 24.0;
01633       }
01634     }
01635     //    return 0.0;
01636     return sum;
01637   }
01638   // Operator applied to a vector U_array
01639   static inline double Local_matrix_bo_cell(const BoCeData* bocedata,
01640                                           double*  U_array, int num_v) {
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 || tet->N1()==num_v  ||
01647          tet->N2()==num_v || tet->N3()==num_v ) {
01648         sum = sum +
01649               (YM0*UAA0 + YM1*UAA1 + YM2*UAA2 + YM3*UAA3) / 24.0;
01650       }
01651     }
01652     //    return 0.0;
01653     return sum;
01654   }
01655   static inline double diag_F_bo_cell(const BoCeData* bocedata,
01656                                       int num_v) {
01657    static double sum;
01658     Tetraeder_storage* tet;
01659 
01660     sum=0.0;
01661     for(tet=bocedata->Give_tets();tet!=NULL;tet=tet->Next()) {
01662       if(tet->N0()==num_v) sum = sum + YM0 / 24.0;
01663       if(tet->N1()==num_v) sum = sum + YM1 / 24.0;
01664       if(tet->N2()==num_v) sum = sum + YM2 / 24.0;
01665       if(tet->N3()==num_v) sum = sum + YM3 / 24.0;
01666     }
01667     return sum;
01668   }
01669 };
01670 
01671 // Formeln fuer helmdy
01672 class helmdy_FE_const {
01673  public:
01674   diff_type typ_v() { return  Dy_type; }; 
01675   diff_type typ_u() { return  Helm_type; };
01676 
01677   helmdy_FE_const() { }
01678   static inline int    ops_interior() { return 12; }
01679   static inline int    ops_diag_interior() { return 0; }
01680   static inline double stencil(double uN,
01681                                double uW, double uM, double uE,
01682                                double uS, double uT, double uD,
01683                                double uND, double uWN, double uWT,
01684                                double uED, double uST, double uES,
01685                                double uEST, double uWND,
01686                                double h_mesh) {
01687     return  -(  uN + uWN + uND + uT + uE 
01688                 - uS - uES - uST - uD - uW)
01689       * h_mesh * h_mesh / 6.0; 
01690   }
01691   static inline double diag_stencil(double h_mesh) {
01692     return 0.0;
01693   }
01694   static inline double F_reg_cell(double*const * u_Recell, int num_var,
01695                                   dir_sons dir_v,
01696                                   double h_mesh, double* sten) {
01697     static int i;
01698     static double sum;
01699   
01700     sum = 0.0;
01701     for(i=0;i<8;++i) sum = sum + sten[dir_v+8*i] * u_Recell[i][num_var];
01702     return sum * h_mesh * h_mesh / 24.0;
01703   }
01704   static inline double Factor_reg_cell(double h_mesh) {
01705     return h_mesh * h_mesh / 24.0;
01706   }
01707   static inline double diag_F_reg_cell(dir_sons dir_v,
01708                                        double h_mesh, double* sten) {
01709     return sten[dir_v*9] * h_mesh * h_mesh / 24.0;
01710   }
01711 
01712   // Operator applied to variable with number
01713   static inline double F_bo_cell(const BoCeData* bocedata,
01714                                  int num_v, int num_variable) {
01715     static double sum;
01716     Tetraeder_storage* tet;
01717 
01718     sum=0.0;
01719     for(tet=bocedata->Give_tets();tet!=NULL;tet=tet->Next()) {
01720       if(tet->N0()==num_v) sum += YM0*(UU0 + UU1 + UU2 + UU3) / 24.0;
01721       if(tet->N1()==num_v) sum += YM1*(UU0 + UU1 + UU2 + UU3) / 24.0;
01722       if(tet->N2()==num_v) sum += YM2*(UU0 + UU1 + UU2 + UU3) / 24.0;
01723       if(tet->N3()==num_v) sum += YM3*(UU0 + UU1 + UU2 + UU3) / 24.0;
01724     }
01725     return sum;
01726   }
01727   // Operator applied to a vector U_array
01728   static inline double Local_matrix_bo_cell(const BoCeData* bocedata,
01729                                           double*  U_array, int num_v) {
01730    static double sum;
01731     Tetraeder_storage* tet;
01732 
01733     sum=0.0;
01734     for(tet=bocedata->Give_tets();tet!=NULL;tet=tet->Next()) {
01735       if(tet->N0()==num_v) sum += YM0*(UAA0 + UAA1 + UAA2 + UAA3) / 24.0;
01736       if(tet->N1()==num_v) sum += YM1*(UAA0 + UAA1 + UAA2 + UAA3) / 24.0;
01737       if(tet->N2()==num_v) sum += YM2*(UAA0 + UAA1 + UAA2 + UAA3) / 24.0;
01738       if(tet->N3()==num_v) sum += YM3*(UAA0 + UAA1 + UAA2 + UAA3) / 24.0;
01739     }
01740     return sum;
01741   }
01742   static inline double diag_F_bo_cell(const BoCeData* bocedata,
01743                                       int num_v) {
01744    static double sum;
01745     Tetraeder_storage* tet;
01746 
01747     sum=0.0;
01748     for(tet=bocedata->Give_tets();tet!=NULL;tet=tet->Next()) {
01749       if(tet->N0()==num_v) sum += YM0 / 24.0;
01750       if(tet->N1()==num_v) sum += YM1 / 24.0;
01751       if(tet->N2()==num_v) sum += YM2 / 24.0;
01752       if(tet->N3()==num_v) sum += YM3 / 24.0;
01753     }
01754     return sum;
01755   }
01756 };
01757 
01758 // Formeln fuer dzhelm
01759 class dzhelm_FE_const {
01760  public:
01761   diff_type typ_u() { return  Dz_type; }; 
01762   diff_type typ_v() { return  Helm_type; };
01763 
01764   dzhelm_FE_const() { }
01765   static inline int    ops_interior() { return 17; }
01766   static inline int    ops_diag_interior() { return 0; }
01767   static inline double stencil(double uN,
01768                                double uW, double uM, double uE,
01769                                double uS, double uT, double uD,
01770                                double uND, double uWN, double uWT,
01771                                double uED, double uST, double uES,
01772                                double uEST, double uWND,
01773                                double h_mesh) {
01774     return  (3.0*(uT - uD) 
01775              + uWT  - uW   + uE  - uED
01776              + uWN  - uWND + uN  - uND
01777              + uEST - uES  + uST - uS)
01778       * h_mesh * h_mesh / 12.0; 
01779   }
01780   static inline double diag_stencil(double h_mesh) {
01781     return 0.0;
01782   }
01783   static inline double F_reg_cell(double*const * u_Recell, int num_var,
01784                                   dir_sons dir_v,
01785                                   double h_mesh, double* sten) {
01786     static int i;
01787     static double sum;
01788   
01789     sum = 0.0;
01790     for(i=0;i<8;++i) sum = sum + sten[dir_v+8*i] * u_Recell[i][num_var];
01791     //  return 0.0;
01792     return sum * h_mesh * h_mesh / 24.0;
01793   }
01794   static inline double Factor_reg_cell(double h_mesh) {
01795     return h_mesh * h_mesh / 24.0;
01796   }
01797   static inline double diag_F_reg_cell(dir_sons dir_v,
01798                                        double h_mesh, double* sten) {
01799     return sten[dir_v*9] * h_mesh * h_mesh / 24.0;
01800   }
01801 
01802   // Operator applied to variable with number
01803   static inline double F_bo_cell(const BoCeData* bocedata,
01804                                  int num_v, int num_variable) {
01805     static double sum;
01806     Tetraeder_storage* tet;
01807 
01808     sum=0.0;
01809     for(tet=bocedata->Give_tets();tet!=NULL;tet=tet->Next()) {
01810       if(tet->N0()==num_v || tet->N1()==num_v  ||
01811          tet->N2()==num_v || tet->N3()==num_v ) {
01812         sum = sum +
01813               (ZM0*UU0 + ZM1*UU1 + ZM2*UU2 + ZM3*UU3) / 24.0;
01814       }
01815     }
01816     // return 0.0;
01817     return sum;
01818   }
01819   // Operator applied to a vector U_array
01820   static inline double Local_matrix_bo_cell(const BoCeData* bocedata,
01821                                           double*  U_array, int num_v) {
01822   static double sum;
01823     Tetraeder_storage* tet;
01824 
01825     sum=0.0;
01826     for(tet=bocedata->Give_tets();tet!=NULL;tet=tet->Next()) {
01827       if(tet->N0()==num_v || tet->N1()==num_v  ||
01828          tet->N2()==num_v || tet->N3()==num_v ) {
01829         sum = sum +
01830               (ZM0*UAA0 + ZM1*UAA1 + ZM2*UAA2 + ZM3*UAA3) / 24.0;
01831       }
01832     }
01833     // return 0.0;
01834     return sum;
01835   }
01836   static inline double diag_F_bo_cell(const BoCeData* bocedata,
01837                                       int num_v) {
01838    static double sum;
01839     Tetraeder_storage* tet;
01840 
01841     sum=0.0;
01842     for(tet=bocedata->Give_tets();tet!=NULL;tet=tet->Next()) {
01843       if(tet->N0()==num_v) sum = sum + ZM0 / 24.0;
01844       if(tet->N1()==num_v) sum = sum + ZM1 / 24.0;
01845       if(tet->N2()==num_v) sum = sum + ZM2 / 24.0;
01846       if(tet->N3()==num_v) sum = sum + ZM3 / 24.0;
01847     }
01848     return sum;
01849   }
01850 };
01851 
01852 // Formeln fuer helmdz
01853 class helmdz_FE_const {
01854  public:
01855   diff_type typ_v() { return  Dz_type; }; 
01856   diff_type typ_u() { return  Helm_type; };
01857 
01858   helmdz_FE_const() { }
01859   static inline int    ops_interior() { return 17; }
01860   static inline int    ops_diag_interior() { return 0; }
01861   static inline double stencil(double uN,
01862                                double uW, double uM, double uE,
01863                                double uS, double uT, double uD,
01864                                double uND, double uWN, double uWT,
01865                                double uED, double uST, double uES,
01866                                double uEST, double uWND,
01867                                double h_mesh) {
01868     return  -(3.0*(uT - uD) 
01869               + uWT  - uW   + uE  - uED
01870               + uWN  - uWND + uN  - uND
01871               + uEST - uES  + uST - uS)
01872       * h_mesh * h_mesh / 12.0; 
01873   }
01874   static inline double diag_stencil(double h_mesh) {
01875     return 0.0;
01876   }
01877   static inline double F_reg_cell(double*const * u_Recell, int num_var,
01878                                   dir_sons dir_v,
01879                                   double h_mesh, double* sten) {
01880     static int i;
01881     static double sum;
01882   
01883     sum = 0.0;
01884     for(i=0;i<8;++i) sum = sum + sten[dir_v+8*i] * u_Recell[i][num_var];
01885     return sum * h_mesh * h_mesh / 24.0;
01886   }
01887   static inline double Factor_reg_cell(double h_mesh) {
01888     return h_mesh * h_mesh / 24.0;
01889   }
01890   static inline double diag_F_reg_cell(dir_sons dir_v,
01891                                        double h_mesh, double* sten) {
01892     return sten[dir_v*9] * h_mesh * h_mesh / 24.0;
01893   }
01894 
01895   // Operator applied to variable with number
01896   static inline double F_bo_cell(const BoCeData* bocedata,
01897                                  int num_v, int num_variable) {
01898     static double sum;
01899     Tetraeder_storage* tet;
01900 
01901     sum=0.0;
01902     for(tet=bocedata->Give_tets();tet!=NULL;tet=tet->Next()) {
01903       if(tet->N0()==num_v) sum += ZM0*(UU0 + UU1 + UU2 + UU3) / 24.0;
01904       if(tet->N1()==num_v) sum += ZM1*(UU0 + UU1 + UU2 + UU3) / 24.0;
01905       if(tet->N2()==num_v) sum += ZM2*(UU0 + UU1 + UU2 + UU3) / 24.0;
01906       if(tet->N3()==num_v) sum += ZM3*(UU0 + UU1 + UU2 + UU3) / 24.0;
01907     }
01908     // return 0.0;
01909     return sum;
01910   }
01911   // Operator applied to a vector U_array
01912   static inline double Local_matrix_bo_cell(const BoCeData* bocedata,
01913                                           double*  U_array, int num_v) {
01914    static double sum;
01915     Tetraeder_storage* tet;
01916 
01917     sum=0.0;
01918     for(tet=bocedata->Give_tets();tet!=NULL;tet=tet->Next()) {
01919       if(tet->N0()==num_v) sum += ZM0*(UAA0 + UAA1 + UAA2 + UAA3) / 24.0;
01920       if(tet->N1()==num_v) sum += ZM1*(UAA0 + UAA1 + UAA2 + UAA3) / 24.0;
01921       if(tet->N2()==num_v) sum += ZM2*(UAA0 + UAA1 + UAA2 + UAA3) / 24.0;
01922       if(tet->N3()==num_v) sum += ZM3*(UAA0 + UAA1 + UAA2 + UAA3) / 24.0;
01923     }
01924     // return 0.0;
01925     return sum;
01926   }
01927   static inline double diag_F_bo_cell(const BoCeData* bocedata,
01928                                       int num_v) {
01929     static double sum;
01930     Tetraeder_storage* tet;
01931 
01932     sum=0.0;
01933     for(tet=bocedata->Give_tets();tet!=NULL;tet=tet->Next()) {
01934       if(tet->N0()==num_v) sum += ZM0 / 24.0;
01935       if(tet->N1()==num_v) sum += ZM1 / 24.0;
01936       if(tet->N2()==num_v) sum += ZM2 / 24.0;
01937       if(tet->N3()==num_v) sum += ZM3 / 24.0;
01938     }
01939     return sum;
01940   }
01941 };
01942 
01943 #endif
01944
Generated on Fri Nov 2 01:25:58 2007 for IPPL by doxygen 1.3.5