src/expde/formulas/loc_sten.cc

Go to the documentation of this file.

//    expde: expression templates for partial differential equations.
//    Copyright (C) 2001  Christoph Pflaum
//    This program is free software; you can redistribute it and/or modify
//    it under the terms of the GNU General Public License as published by
//    the Free Software Foundation; either version 2 of the License, or
//    (at your option) any later version.
//
//    This program is distributed in the hope that it will be useful,
//    but WITHOUT ANY WARRANTY; without even the implied warranty of
//    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
//    GNU General Public License for more details.
//
//                 SEE  Notice1.doc made by 
//                 LAWRENCE LIVERMORE NATIONAL LABORATORY
//

// ------------------------------------------------------------
//
// loc_sten.cc
//
// ------------------------------------------------------------

// Include level 0:
#ifdef COMP_GNUOLD
 #include <iostream.h>
 #include <fstream.h>
 #include <math.h>
#else
 #include <iostream>
 #include <fstream>
 #include <cmath>
#endif 

// Include level 1:
#include "../paramete.h"
#include "../abbrevi.h"
#include "../math_lib/math_lib.h"

// Include level special
#include "../formulas/loc_sten.h"

#define x0 V0.x
#define y0 V0.y
#define z0 V0.z

#define x1 V1.x
#define y1 V1.y
#define z1 V1.z

#define x2 V2.x
#define y2 V2.y
#define z2 V2.z

#define x3 V3.x
#define y3 V3.y
#define z3 V3.z

inline bool In_tet(dir_sons* tet, dir_sons s) {
  return (tet[0]==s || tet[1]==s || tet[2]==s || tet[3]==s);
}

inline double Tetraeder_Det_dx0(const D3vector& V0, const D3vector& V1,
                                const D3vector& V2, const D3vector& V3) {
  return y3 * (z1 - z2) + y1 * (z2 - z3) + y2 * (-z1 + z3);
}

inline double Tetraeder_Det_dy0(const D3vector& V0, const D3vector& V1,
                                const D3vector& V2, const D3vector& V3) {
  return x3 * (-z1 + z2) + x2 * (z1 - z3) + x1 * (-z2 + z3);
}

inline double Tetraeder_Det_dz0(const D3vector& V0, const D3vector& V1,
                                const D3vector& V2, const D3vector& V3) {
  return x3 * (y1 - y2) + x1 * (y2 - y3) + x2 * (-y1 + y3);
}

inline double Tetraeder_Det_dx(int n, 
                               const dir_sons co0, const dir_sons co1,
                               const dir_sons co2, const dir_sons co3) {
  D3vector V0, V1, V2, V3;
  x0 = x1DCoord(co0);  y0 = y1DCoord(co0);  z0 = z1DCoord(co0);
  x1 = x1DCoord(co1);  y1 = y1DCoord(co1);  z1 = z1DCoord(co1);
  x2 = x1DCoord(co2);  y2 = y1DCoord(co2);  z2 = z1DCoord(co2);
  x3 = x1DCoord(co3);  y3 = y1DCoord(co3);  z3 = z1DCoord(co3);

  if(n == 0) return Tetraeder_Det_dx0(V0,V1,V2,V3);
  if(n == 1) return Tetraeder_Det_dx0(V1,V3,V2,V0);
  if(n == 2) return Tetraeder_Det_dx0(V2,V1,V3,V0);
  //  if(n == 3) 
  return Tetraeder_Det_dx0(V3,V0,V2,V1);
};

inline double Tetraeder_Det_dy(int n, 
                               const dir_sons co0, const dir_sons co1,
                               const dir_sons co2, const dir_sons co3) {
  D3vector V0, V1, V2, V3;
  x0 = x1DCoord(co0);  y0 = y1DCoord(co0);  z0 = z1DCoord(co0);
  x1 = x1DCoord(co1);  y1 = y1DCoord(co1);  z1 = z1DCoord(co1);
  x2 = x1DCoord(co2);  y2 = y1DCoord(co2);  z2 = z1DCoord(co2);
  x3 = x1DCoord(co3);  y3 = y1DCoord(co3);  z3 = z1DCoord(co3);

  if(n == 0) return Tetraeder_Det_dy0(V0,V1,V2,V3);
  if(n == 1) return Tetraeder_Det_dy0(V1,V3,V2,V0);
  if(n == 2) return Tetraeder_Det_dy0(V2,V1,V3,V0);
  //  if(n == 3) 
  return Tetraeder_Det_dy0(V3,V0,V2,V1);
};

inline double Tetraeder_Det_dz(int n, 
                               const dir_sons co0, const dir_sons co1,
                               const dir_sons co2, const dir_sons co3) {
  D3vector V0, V1, V2, V3;
  x0 = x1DCoord(co0);  y0 = y1DCoord(co0);  z0 = z1DCoord(co0);
  x1 = x1DCoord(co1);  y1 = y1DCoord(co1);  z1 = z1DCoord(co1);
  x2 = x1DCoord(co2);  y2 = y1DCoord(co2);  z2 = z1DCoord(co2);
  x3 = x1DCoord(co3);  y3 = y1DCoord(co3);  z3 = z1DCoord(co3);

  if(n == 0) return Tetraeder_Det_dz0(V0,V1,V2,V3);
  if(n == 1) return Tetraeder_Det_dz0(V1,V3,V2,V0);
  if(n == 2) return Tetraeder_Det_dz0(V2,V1,V3,V0);
  //  if(n == 3) 
  return Tetraeder_Det_dz0(V3,V0,V2,V1);
};


inline double Evaluate_tet(dir_sons* tet, dir_sons s, diff_type dt) {
  if(dt==Dx_type) {
    if(tet[0]==s) return Tetraeder_Det_dx(0,tet[0],tet[1],tet[2],tet[3]);
    if(tet[1]==s) return Tetraeder_Det_dx(1,tet[0],tet[1],tet[2],tet[3]);
    if(tet[2]==s) return Tetraeder_Det_dx(2,tet[0],tet[1],tet[2],tet[3]);
    // if(tet[3]==s) 
    return Tetraeder_Det_dx(3,tet[0],tet[1],tet[2],tet[3]);
  }
  if(dt==Dy_type) {
    if(tet[0]==s) return Tetraeder_Det_dy(0,tet[0],tet[1],tet[2],tet[3]);
    if(tet[1]==s) return Tetraeder_Det_dy(1,tet[0],tet[1],tet[2],tet[3]);
    if(tet[2]==s) return Tetraeder_Det_dy(2,tet[0],tet[1],tet[2],tet[3]);
    // if(tet[3]==s) 
    return Tetraeder_Det_dy(3,tet[0],tet[1],tet[2],tet[3]);
  }
  //  if(dt==Dz_type) {
  if(tet[0]==s) return Tetraeder_Det_dz(0,tet[0],tet[1],tet[2],tet[3]);
  if(tet[1]==s) return Tetraeder_Det_dz(1,tet[0],tet[1],tet[2],tet[3]);
  if(tet[2]==s) return Tetraeder_Det_dz(2,tet[0],tet[1],tet[2],tet[3]);
  // if(tet[3]==s) 
  return Tetraeder_Det_dz(3,tet[0],tet[1],tet[2],tet[3]);
  //  }
}



void Calc_Loc_stencil(double* sten, diff_type typ_u, diff_type typ_v) {
  int u,v, k;
  dir_sons Tets[6][4];
  
  /* Alternative Konstruktion:
  // 0. Tet: WNT, END, WST, EST - schraeg
  Tets[0][0] = WNTd;
  Tets[0][1] = ENDd;
  Tets[0][2] = WSTd; 
  Tets[0][3] = ESTd;     

  // 1. Tet: EST, END, WST, ESD - Seite    
  Tets[1][0] = ESTd;
  Tets[1][1] = ENDd;
  Tets[1][2] = WSTd; 
  Tets[1][3] = ESDd;     

  // 2. Tet: WND, WSD, WST, ESD - Ecke  
  Tets[2][0] = WNDd;
  Tets[2][1] = WSDd;
  Tets[2][2] = WSTd; 
  Tets[2][3] = ESDd;     

  // 3. Tet: WST, WND, ESD, END -
  Tets[3][0] = WSTd;
  Tets[3][1] = WNDd;
  Tets[3][2] = ESDd; 
  Tets[3][3] = ENDd;     

  // 4. Tet: ENT, WNT, EST, END - Ecke
  Tets[4][0] = ENTd;
  Tets[4][1] = WNTd;
  Tets[4][2] = ESTd; 
  Tets[4][3] = ENDd;     

  // 5. Tet: WNT, WND, WST, END - Seite
  Tets[5][0] = WNTd;
  Tets[5][1] = WNDd;
  Tets[5][2] = WSTd; 
  Tets[5][3] = ENDd;     
  */
  
  // 0. Tet: WNT, WND, WST, EST 
  Tets[0][0] = WNTd;
  Tets[0][1] = WNDd;
  Tets[0][2] = WSTd; 
  Tets[0][3] = ESTd;     

  // 1. Tet: EST, WND, WST, ESD     
  Tets[1][0] = ESTd;
  Tets[1][1] = WNDd;
  Tets[1][2] = WSTd; 
  Tets[1][3] = ESDd;     

  // 2. Tet: WND, WSD, WST, ESD        
  Tets[2][0] = WNDd;
  Tets[2][1] = WSDd;
  Tets[2][2] = WSTd; 
  Tets[2][3] = ESDd;     

  // 3. Tet: EST, WND, ESD, END
  Tets[3][0] = ESTd;
  Tets[3][1] = WNDd;
  Tets[3][2] = ESDd; 
  Tets[3][3] = ENDd;     

  // 4. Tet: ENT, WNT, EST, END
  Tets[4][0] = ENTd;
  Tets[4][1] = WNTd;
  Tets[4][2] = ESTd; 
  Tets[4][3] = ENDd;     

  // 5. Tet: WNT, WND, EST, END
  Tets[5][0] = WNTd;
  Tets[5][1] = WNDd;
  Tets[5][2] = ESTd; 
  Tets[5][3] = ENDd;     
  

  for(u=0;u<8;++u)
    for(v=0;v<8;++v) sten[v+8*u] = 0.0;

  if(typ_u < Helm_type && typ_v < Helm_type)
    for(v=0;v<8;++v)
      for(u=0;u<8;++u) {
        for(k=0;k<6;++k) {
          if(In_tet(Tets[k],(dir_sons)u) && In_tet(Tets[k],(dir_sons)v)) {
            sten[v+8*u] = sten[v+8*u] +
              Evaluate_tet(Tets[k],(dir_sons)u,typ_u) *
              Evaluate_tet(Tets[k],(dir_sons)v,typ_v);
          }
        }
      }

  if(typ_u == Helm_type && typ_v == Helm_type)
    for(v=0;v<8;++v)
      for(u=0;u<8;++u) {
        for(k=0;k<6;++k) {
          if(In_tet(Tets[k],(dir_sons)u) && In_tet(Tets[k],(dir_sons)v)) {
            if(v==u) sten[v+8*u] = sten[v+8*u] + 2;
            else     sten[v+8*u] = sten[v+8*u] + 1;
          }
        }
      }

  if(typ_u < Helm_type && typ_v == Helm_type)
    for(v=0;v<8;++v)
      for(u=0;u<8;++u) {
        for(k=0;k<6;++k) {
          if(In_tet(Tets[k],(dir_sons)u) && In_tet(Tets[k],(dir_sons)v)) {
            sten[v+8*u] = sten[v+8*u] +
              Evaluate_tet(Tets[k],(dir_sons)u,typ_u);
          }
        }
      }
  if(typ_u == Helm_type && typ_v < Helm_type)
    for(v=0;v<8;++v)
      for(u=0;u<8;++u) {
        for(k=0;k<6;++k) {
          if(In_tet(Tets[k],(dir_sons)u) && In_tet(Tets[k],(dir_sons)v)) {
            sten[v+8*u] = sten[v+8*u] +
              Evaluate_tet(Tets[k],(dir_sons)v,typ_v);
          }
        }
      }
}

void Print_u_Recell(double*const * u_Recell, int num) {
  cout << "\n  North: " << endl;
  cout << "    " << u_Recell[WNTd][num] << ", " << u_Recell[ENTd][num] << endl;
  cout << "    " << u_Recell[WNDd][num] << ", " << u_Recell[ENDd][num] << endl;
  cout << "  South: " << endl;
  cout << "    " << u_Recell[WSTd][num] << ", " << u_Recell[ESTd][num] << endl;
  cout << "    " << u_Recell[WSDd][num] << ", " << u_Recell[ESDd][num] << endl;
}

void Print_part_loc_stencil(double* sten) {
  cout << "\n   North: " << endl;
  cout << "     " << sten[WNTd] << ", " << sten[ENTd] << endl;
  cout << "     " << sten[WNDd] << ", " << sten[ENDd] << endl;
  cout << "   South: " << endl;
  cout << "     " << sten[WSTd] << ", " << sten[ESTd] << endl;
  cout << "     " << sten[WSDd] << ", " << sten[ESDd] << endl;
}

void Print_Loc_stencil(double* sten) {
  int i;
  for(i=0;i<8;++i) {
    cout << "\n Sten bei: ";  Print((dir_sons)i);
    Print_part_loc_stencil(&(sten[8*i]));
  }
}

Loc_stencil_container::Loc_stencil_container() {
  int i, j;
  sten = new double**[4];
  for(i=0;i<4;++i) sten[i] = new double*[4];  
  for(i=0;i<4;++i) for(j=0;j<4;++j) sten[i][j] = NULL;
  sten_Poisson = new double[64];

  // Poisson berechnen
  sten_sp = Give_Loc_sten(Dx_type,Dx_type);
  for(i=0;i<8;++i) for(j=0;j<8;++j)
    sten_Poisson[i+8*j] = sten_sp[i+8*j];
  sten_sp = Give_Loc_sten(Dy_type,Dy_type);
  for(i=0;i<8;++i) for(j=0;j<8;++j)
    sten_Poisson[i+8*j] = sten_Poisson[i+8*j] + sten_sp[i+8*j];
  sten_sp = Give_Loc_sten(Dz_type,Dz_type);
  for(i=0;i<8;++i) for(j=0;j<8;++j)
    sten_Poisson[i+8*j] = sten_Poisson[i+8*j] + sten_sp[i+8*j];
};

double* Loc_stencil_container::Give_Loc_sten(diff_type typ_u,
                                              diff_type typ_v) const {
  if(sten[typ_u][typ_v] == NULL) {
    sten[typ_u][typ_v] = new double[64];
    Calc_Loc_stencil(sten[typ_u][typ_v],typ_u,typ_v);
  }
  return sten[typ_u][typ_v];
};

void Print_27_stencil(diff_type typ_u, diff_type typ_v) {
  int i,j,k, u,v;
  double Sten27[27];
  double* sten;
  sten = new double[64];
  for(i=0;i<27;++i) Sten27[i] = 0;

  Calc_Loc_stencil(sten, typ_u, typ_v);
  
  Print_Factor(typ_u,typ_v);

  // Nachbarzelle: j
  for(i=0;i<8;++i) for(j=0;j<8;++j) {
    v = opposite3D((dir_sons)j);
    u = i;
    Sten27[trans_cellP_27((dir_sons)i,(dir_sons)j)] = 
      Sten27[trans_cellP_27((dir_sons)i,(dir_sons)j)] + sten[v+8*u];
  }
  j=2;
  cout << "\n Nord: \n";
  for(k=2;k>=0;--k) {
    for(i=0;i<3;++i) {
      cout << ",   " << Sten27[trans27((Ort1D)i,(Ort1D)j,(Ort1D)k)];
    }
    cout << "\n";
  }
  j=1;
  cout << "\n Mitte: \n";
  for(k=2;k>=0;--k) {
    for(i=0;i<3;++i) {
      cout << ",   " << Sten27[trans27((Ort1D)i,(Ort1D)j,(Ort1D)k)];
    }
    cout << "\n";
  }
  j=0;
  cout << "\n Sued: \n";
  for(k=2;k>=0;--k) {
    for(i=0;i<3;++i) {
      cout << ",   " << Sten27[trans27((Ort1D)i,(Ort1D)j,(Ort1D)k)];
    }
    cout << "\n";
  }
}

void Print_27_stencil_var(diff_type typ_u, diff_type typ_v) {
  int i,j,k, u,v, s;
  double Sten27[27];
  double* sten;
  sten = new double[64];
  for(i=0;i<27;++i) Sten27[i] = 0;

  Calc_Loc_stencil(sten, typ_u, typ_v);
  
  Print_Factor(typ_u,typ_v);

  // Nachbarzelle: s
  for(s=0;s<8;++s) {
    cout << " cell: "; Print((dir_sons)s);  cout << endl;
    for(i=0;i<27;++i) Sten27[i] = 0;
    for(i=0;i<8;++i) {
      v = opposite3D((dir_sons)s);
      u = i;
      Sten27[trans_cellP_27((dir_sons)i,(dir_sons)s)] = 
        Sten27[trans_cellP_27((dir_sons)i,(dir_sons)s)] + sten[v+8*u];
    }
    j=2;
    cout << "\n Nord: \n";
    for(k=2;k>=0;--k) {
      for(i=0;i<3;++i) {
        cout << ",   " << Sten27[trans27((Ort1D)i,(Ort1D)j,(Ort1D)k)];
      }
      cout << "\n";
    }
    j=1;
    cout << "\n Mitte: \n";
    for(k=2;k>=0;--k) {
      for(i=0;i<3;++i) {
        cout << ",   " << Sten27[trans27((Ort1D)i,(Ort1D)j,(Ort1D)k)];
      }
      cout << "\n";
    }
    j=0;
    cout << "\n Sued: \n";
    for(k=2;k>=0;--k) {
      for(i=0;i<3;++i) {
        cout << ",   " << Sten27[trans27((Ort1D)i,(Ort1D)j,(Ort1D)k)];
      }
      cout << "\n";
    }
  }
}



void Print_with_Vorzeichen(double x) {
  if(ABS(x-1) < 0.001) {
    cout << " + ";
  }
  else if(ABS(x+1) < 0.001) {
    cout << " - ";
  }
  else if(x>0.1) cout << " + " << x << "*";
  else cout << " - " << -x << "*";
}



void Print_27_stencil_formula_var(diff_type typ_u, diff_type typ_v) {
  int i,j,k, u,v, s;
  double Sten27[27];
  double* sten;
  
  cout.precision(1);

  sten = new double[64];
  for(i=0;i<27;++i) Sten27[i] = 0;

  Calc_Loc_stencil(sten, typ_u, typ_v);
  
  Print_Factor(typ_u,typ_v);

  // Nachbarzelle: s
  for(s=0;s<8;++s) {
    //    cout << " cell: "; Print((dir_sons)s);
    //    cout << " ------------------------------------ ";
    cout << " \n var_interpo_"; Print((dir_sons)s); cout << " * (";

    for(i=0;i<27;++i) Sten27[i] = 0;
    for(i=0;i<8;++i) {
      v = opposite3D((dir_sons)s);
      u = i;
      Sten27[trans_cellP_27((dir_sons)i,(dir_sons)s)] = 
        Sten27[trans_cellP_27((dir_sons)i,(dir_sons)s)] + sten[v+8*u];
    }
    // Nord
    // ----
    j=2;

    i=0;  // West
    k=2; 
    if(ABS(Sten27[trans27((Ort1D)i,(Ort1D)j,(Ort1D)k)]) > 0.01) {
      Print_with_Vorzeichen(Sten27[trans27((Ort1D)i,(Ort1D)j,(Ort1D)k)]);;
      cout << "uWNT";
    }
    k=1;
    if(ABS(Sten27[trans27((Ort1D)i,(Ort1D)j,(Ort1D)k)]) > 0.01) {
      Print_with_Vorzeichen(Sten27[trans27((Ort1D)i,(Ort1D)j,(Ort1D)k)]);;
      cout << "uWN";
    }
    k=0;
    if(ABS(Sten27[trans27((Ort1D)i,(Ort1D)j,(Ort1D)k)]) > 0.01) {
      Print_with_Vorzeichen(Sten27[trans27((Ort1D)i,(Ort1D)j,(Ort1D)k)]);;
      cout << "uWND";
    }

    i=1;  // Middle
    k=2; 
    if(ABS(Sten27[trans27((Ort1D)i,(Ort1D)j,(Ort1D)k)]) > 0.01) {
      Print_with_Vorzeichen(Sten27[trans27((Ort1D)i,(Ort1D)j,(Ort1D)k)]);;
      cout << "uNT";
    }
    k=1;
    if(ABS(Sten27[trans27((Ort1D)i,(Ort1D)j,(Ort1D)k)]) > 0.01) {
      Print_with_Vorzeichen(Sten27[trans27((Ort1D)i,(Ort1D)j,(Ort1D)k)]);;
      cout << "uN";
    }
    k=0;
    if(ABS(Sten27[trans27((Ort1D)i,(Ort1D)j,(Ort1D)k)]) > 0.01) {
      Print_with_Vorzeichen(Sten27[trans27((Ort1D)i,(Ort1D)j,(Ort1D)k)]);;
      cout << "uND";
    }

    i=2;  // East
    k=2; 
    if(ABS(Sten27[trans27((Ort1D)i,(Ort1D)j,(Ort1D)k)]) > 0.01) {
      Print_with_Vorzeichen(Sten27[trans27((Ort1D)i,(Ort1D)j,(Ort1D)k)]);;
      cout << "uENT";
    }
    k=1;
    if(ABS(Sten27[trans27((Ort1D)i,(Ort1D)j,(Ort1D)k)]) > 0.01) {
      Print_with_Vorzeichen(Sten27[trans27((Ort1D)i,(Ort1D)j,(Ort1D)k)]);;
      cout << "uEN";
    }
    k=0;
    if(ABS(Sten27[trans27((Ort1D)i,(Ort1D)j,(Ort1D)k)]) > 0.01) {
      Print_with_Vorzeichen(Sten27[trans27((Ort1D)i,(Ort1D)j,(Ort1D)k)]);;
      cout << "uEND";
    }

    // Middle
    // ------
    j=1;

    i=0;  // West
    k=2; 
    if(ABS(Sten27[trans27((Ort1D)i,(Ort1D)j,(Ort1D)k)]) > 0.01) {
      Print_with_Vorzeichen(Sten27[trans27((Ort1D)i,(Ort1D)j,(Ort1D)k)]);;
      cout << "uWT";
    }
    k=1;
    if(ABS(Sten27[trans27((Ort1D)i,(Ort1D)j,(Ort1D)k)]) > 0.01) {
      Print_with_Vorzeichen(Sten27[trans27((Ort1D)i,(Ort1D)j,(Ort1D)k)]);;
      cout << "uW";
    }
    k=0;
    if(ABS(Sten27[trans27((Ort1D)i,(Ort1D)j,(Ort1D)k)]) > 0.01) {
      Print_with_Vorzeichen(Sten27[trans27((Ort1D)i,(Ort1D)j,(Ort1D)k)]);;
      cout << "uWD";
    }

    i=1;  // Middle
    k=2; 
    if(ABS(Sten27[trans27((Ort1D)i,(Ort1D)j,(Ort1D)k)]) > 0.01) {
      Print_with_Vorzeichen(Sten27[trans27((Ort1D)i,(Ort1D)j,(Ort1D)k)]);;
      cout << "uT";
    }
    k=1;
    if(ABS(Sten27[trans27((Ort1D)i,(Ort1D)j,(Ort1D)k)]) > 0.01) {
      Print_with_Vorzeichen(Sten27[trans27((Ort1D)i,(Ort1D)j,(Ort1D)k)]);;
      cout << "uM";
    }
    k=0;
    if(ABS(Sten27[trans27((Ort1D)i,(Ort1D)j,(Ort1D)k)]) > 0.01) {
      Print_with_Vorzeichen(Sten27[trans27((Ort1D)i,(Ort1D)j,(Ort1D)k)]);;
      cout << "uD";
    }

    i=2;  // East
    k=2; 
    if(ABS(Sten27[trans27((Ort1D)i,(Ort1D)j,(Ort1D)k)]) > 0.01) {
      Print_with_Vorzeichen(Sten27[trans27((Ort1D)i,(Ort1D)j,(Ort1D)k)]);;
      cout << "uET";
    }
    k=1;
    if(ABS(Sten27[trans27((Ort1D)i,(Ort1D)j,(Ort1D)k)]) > 0.01) {
      Print_with_Vorzeichen(Sten27[trans27((Ort1D)i,(Ort1D)j,(Ort1D)k)]);;
      cout << "uE";
    }
    k=0;
    if(ABS(Sten27[trans27((Ort1D)i,(Ort1D)j,(Ort1D)k)]) > 0.01) {
      Print_with_Vorzeichen(Sten27[trans27((Ort1D)i,(Ort1D)j,(Ort1D)k)]);;
      cout << "uED";
    }

    //Sued
    // ----
    j=0;

    i=0;  // West
    k=2; 
    if(ABS(Sten27[trans27((Ort1D)i,(Ort1D)j,(Ort1D)k)]) > 0.01) {
      Print_with_Vorzeichen(Sten27[trans27((Ort1D)i,(Ort1D)j,(Ort1D)k)]);;
      cout << "uWST";
    }
    k=1;
    if(ABS(Sten27[trans27((Ort1D)i,(Ort1D)j,(Ort1D)k)]) > 0.01) {
      Print_with_Vorzeichen(Sten27[trans27((Ort1D)i,(Ort1D)j,(Ort1D)k)]);;
      cout << "uWS";
    }
    k=0;
    if(ABS(Sten27[trans27((Ort1D)i,(Ort1D)j,(Ort1D)k)]) > 0.01) {
      Print_with_Vorzeichen(Sten27[trans27((Ort1D)i,(Ort1D)j,(Ort1D)k)]);;
      cout << "uWSD";
    }

    i=1;  // Middle
    k=2; 
    if(ABS(Sten27[trans27((Ort1D)i,(Ort1D)j,(Ort1D)k)]) > 0.01) {
      Print_with_Vorzeichen(Sten27[trans27((Ort1D)i,(Ort1D)j,(Ort1D)k)]);;
      cout << "uST";
    }
    k=1;
    if(ABS(Sten27[trans27((Ort1D)i,(Ort1D)j,(Ort1D)k)]) > 0.01) {
      Print_with_Vorzeichen(Sten27[trans27((Ort1D)i,(Ort1D)j,(Ort1D)k)]);;
      cout << "uS";
    }
    k=0;
    if(ABS(Sten27[trans27((Ort1D)i,(Ort1D)j,(Ort1D)k)]) > 0.01) {
      Print_with_Vorzeichen(Sten27[trans27((Ort1D)i,(Ort1D)j,(Ort1D)k)]);;
      cout << "uSD";
    }

    i=2;  // East
    k=2; 
    if(ABS(Sten27[trans27((Ort1D)i,(Ort1D)j,(Ort1D)k)]) > 0.01) {
      Print_with_Vorzeichen(Sten27[trans27((Ort1D)i,(Ort1D)j,(Ort1D)k)]);;
      cout << "uEST";
    }
    k=1;
    if(ABS(Sten27[trans27((Ort1D)i,(Ort1D)j,(Ort1D)k)]) > 0.01) {
      Print_with_Vorzeichen(Sten27[trans27((Ort1D)i,(Ort1D)j,(Ort1D)k)]);;
      cout << "uES";
    }
    k=0;
    if(ABS(Sten27[trans27((Ort1D)i,(Ort1D)j,(Ort1D)k)]) > 0.01) {
      Print_with_Vorzeichen(Sten27[trans27((Ort1D)i,(Ort1D)j,(Ort1D)k)]);;
      cout << "uESD";
    }
    cout << ") + ";
  }
  cout << "\n Formula of stencil!" << endl;
}




void Print_27_stencil_of_local_matrix(double* sten) {
  int i,j,k, u,v;
  double Sten27[27];

  for(i=0;i<27;++i) Sten27[i] = 0;

  // Nachbarzelle: j
  for(i=0;i<8;++i) for(j=0;j<8;++j) {
    v = opposite3D((dir_sons)j);
    u = i;
    Sten27[trans_cellP_27((dir_sons)i,(dir_sons)j)] = 
      Sten27[trans_cellP_27((dir_sons)i,(dir_sons)j)] + sten[v+8*u];
  }
  j=2;
  cout << "\n Nord: \n";
  for(k=2;k>=0;--k) {
    for(i=0;i<3;++i) {
      cout << ",   " << Sten27[trans27((Ort1D)i,(Ort1D)j,(Ort1D)k)];
    }
    cout << "\n";
  }
  j=1;
  cout << "\n Mitte: \n";
  for(k=2;k>=0;--k) {
    for(i=0;i<3;++i) {
      cout << ",   " << Sten27[trans27((Ort1D)i,(Ort1D)j,(Ort1D)k)];
    }
    cout << "\n";
  }
  j=0;
  cout << "\n Sued: \n";
  for(k=2;k>=0;--k) {
    for(i=0;i<3;++i) {
      cout << ",   " << Sten27[trans27((Ort1D)i,(Ort1D)j,(Ort1D)k)];
    }
    cout << "\n";
  }
}

void Print_27_stencil(double Sten27[27]) {
  int i,j,k;

  j=2;
  cout << "\n Nord: \n";
  for(k=2;k>=0;--k) {
    for(i=0;i<3;++i) {
      cout << ",   " << Sten27[trans27((Ort1D)i,(Ort1D)j,(Ort1D)k)];
    }
    cout << "\n";
  }
  j=1;
  cout << "\n Mitte: \n";
  for(k=2;k>=0;--k) {
    for(i=0;i<3;++i) {
      cout << ",   " << Sten27[trans27((Ort1D)i,(Ort1D)j,(Ort1D)k)];
    }
    cout << "\n";
  }
  j=0;
  cout << "\n Sued: \n";
  for(k=2;k>=0;--k) {
    for(i=0;i<3;++i) {
      cout << ",   " << Sten27[trans27((Ort1D)i,(Ort1D)j,(Ort1D)k)];
    }
    cout << "\n";
  }
}


void Print_27_stencil_Poisson() {
  int i,j,k, u,v;
  double Sten27[27];
  double sten[64];
  double sten_sp[64];
  /*
  sten = new double*[8];
  for(i=0;i<8;++i) sten[i] = new double[8];
  for(i=0;i<27;++i) Sten27[i] = 0;
  sten_sp = new double*[8];
  for(i=0;i<8;++i) sten_sp[i] = new double[8];
  */
  for(i=0;i<27;++i) Sten27[i] = 0;

  Calc_Loc_stencil(sten, Dx_type, Dx_type);
  for(i=0;i<8;++i) for(j=0;j<8;++j)
  Calc_Loc_stencil(sten_sp, Dy_type, Dy_type);
  for(i=0;i<8;++i) for(j=0;j<8;++j)
    sten[i+8*j] = sten[i+8*j] + sten_sp[i+8*j];
  Calc_Loc_stencil(sten_sp, Dz_type, Dz_type);
  for(i=0;i<8;++i) for(j=0;j<8;++j)
    sten[i+8*j] = sten[i+8*j] + sten_sp[i+8*j];

  Print_Factor_Poisson();

  // Nachbarzelle: j
  for(i=0;i<8;++i) for(j=0;j<8;++j) {
    v = opposite3D((dir_sons)j);
    u = i;
    Sten27[trans_cellP_27((dir_sons)i,(dir_sons)j)] = 
      Sten27[trans_cellP_27((dir_sons)i,(dir_sons)j)] + sten[v+8*u];
  }
  j=2;
  cout << "\n Nord: \n";
  for(k=2;k>=0;--k) {
    for(i=0;i<3;++i) {
      cout << ",   " << Sten27[trans27((Ort1D)i,(Ort1D)j,(Ort1D)k)];
    }
    cout << "\n";
  }
  j=1;
  cout << "\n Mitte: \n";
  for(k=2;k>=0;--k) {
    for(i=0;i<3;++i) {
      cout << ",   " << Sten27[trans27((Ort1D)i,(Ort1D)j,(Ort1D)k)];
    }
    cout << "\n";
  }
  j=0;
  cout << "\n Sued: \n";
  for(k=2;k>=0;--k) {
    for(i=0;i<3;++i) {
      cout << ",   " << Sten27[trans27((Ort1D)i,(Ort1D)j,(Ort1D)k)];
    }
    cout << "\n";
  }
}


void Print_27_stencil_Poisson_var() {
  int i,j,k, u,v,s;
  double Sten27[27];
  double sten[64];
  double sten_sp[64];

  for(i=0;i<27;++i) Sten27[i] = 0;

  Calc_Loc_stencil(sten, Dx_type, Dx_type);
  for(i=0;i<8;++i) for(j=0;j<8;++j)
  Calc_Loc_stencil(sten_sp, Dy_type, Dy_type);
  for(i=0;i<8;++i) for(j=0;j<8;++j)
    sten[i+8*j] = sten[i+8*j] + sten_sp[i+8*j];
  Calc_Loc_stencil(sten_sp, Dz_type, Dz_type);
  for(i=0;i<8;++i) for(j=0;j<8;++j)
    sten[i+8*j] = sten[i+8*j] + sten_sp[i+8*j];

  Print_Factor_Poisson();

  // Nachbarzelle: s
  for(s=0;s<8;++s) {
    cout << " cell: "; Print((dir_sons)s);  cout << endl;
    for(i=0;i<27;++i) Sten27[i] = 0;
    for(i=0;i<8;++i) {
      v = opposite3D((dir_sons)s);
      u = i;
      Sten27[trans_cellP_27((dir_sons)i,(dir_sons)s)] = 
        Sten27[trans_cellP_27((dir_sons)i,(dir_sons)s)] + sten[v+8*u];
    }

    j=2;
    cout << "\n Nord: \n";
    for(k=2;k>=0;--k) {
      for(i=0;i<3;++i) {
        cout << ",   " << Sten27[trans27((Ort1D)i,(Ort1D)j,(Ort1D)k)];
      }
      cout << "\n";
    }
    j=1;
    cout << "\n Mitte: \n";
    for(k=2;k>=0;--k) {
      for(i=0;i<3;++i) {
        cout << ",   " << Sten27[trans27((Ort1D)i,(Ort1D)j,(Ort1D)k)];
      }
      cout << "\n";
    }
    j=0;
    cout << "\n Sued: \n";
    for(k=2;k>=0;--k) {
      for(i=0;i<3;++i) {
        cout << ",   " << Sten27[trans27((Ort1D)i,(Ort1D)j,(Ort1D)k)];
      }
      cout << "\n";
    }
  }
}



void Print_Factor(diff_type typ_u, diff_type typ_v) {
  if(typ_u < Helm_type && typ_v < Helm_type) {
    cout << "\n Faktor ist: h_mesh / 6.0; " << endl;
    return;
  }
  if(typ_u == Helm_type && typ_v == Helm_type) {
    cout << "\n Faktor ist: h_mesh^3 / 120.0; " << endl;
    return;
  }
  cout << "\n Faktor ist: h_mesh^2 / 24.0; " << endl;
}

void Print_Factor_Poisson() {
  cout << "\n Faktor ist: h_mesh / 6.0; " << endl;
}

double Test_evaluation_of_lsm(double* sten, dir_sons dir_v,
                              bool* corner_in_domain) {
  static int i;
  static double sum;
  sum = 0.0;
  for(i=0;i<8;++i) {
    if(corner_in_domain[i])
      sum = sum + sten[dir_v+8*i] * 1.0;
  }
  return sum;
}

double Test_evaluation_of_lsm(double* sten, dir_sons dir_v,
                              double* UU) {
  static int i;
  static double sum;
  sum = 0.0;
  for(i=0;i<8;++i) {
    sum = sum + sten[dir_v+8*i] * UU[i];
  }
  return sum;
}

double Test_evaluation_of_lsm(double* sten, dir_sons dir_v) {
  static int i;
  static double sum;
  sum = 0.0;
  for(i=0;i<8;++i) {
    sum = sum + sten[dir_v+8*i] * 1.0;
  }
  return sum;
}


// for restriction operator

Restriction_stencil_container::Restriction_stencil_container() {
  int i,j,k;

  for(i=0;i<5;i++) for(j=0;j<5;j++) for(k=0;k<5;k++) 
    reg_formula[i][j][k] = 0;
  
  reg_formula[2][2][2] = 2;
  reg_formula[2][3][2] = 1;
  reg_formula[2][1][2] = 1;
  reg_formula[3][2][2] = 1;
  reg_formula[1][2][2] = 1;
  reg_formula[2][2][3] = 1;
  reg_formula[2][2][1] = 1;

  reg_formula[2][3][1] = 1;
  reg_formula[1][3][2] = 1;
  reg_formula[1][2][3] = 1;
  reg_formula[3][2][1] = 1;
  reg_formula[2][1][3] = 1;
  reg_formula[3][1][2] = 1;

  reg_formula[3][1][3] = 1;
  reg_formula[1][3][1] = 1;
}

// #define vom_feinsten_Gitter false
#define vom_feinsten_Gitter true
#define Drucke_code true

void Print_ijk_name(int i, int j, int k) {
  if(Drucke_code==false) {
    cout << " ";
    if(i==2 && j==2 && k==2) 
      cout << "M";
    else {
      if(i==0) cout << "WW";
      if(i==1) cout << "W";
      if(i==3) cout << "E";
      if(i==4) cout << "EE";
      
      if(j==0) cout << "SS";
      if(j==1) cout << "S";
      if(j==3) cout << "N";
      if(j==4) cout << "NN";
      
      if(k==0) cout << "DD";
      if(k==1) cout << "D";
      if(k==3) cout << "T";
      if(k==4) cout << "TT";
    }
    cout << " ";
  }
  else {
    if(i==2 && j==2 && k==2) {}
    else {
      if(i==0) cout << ".next_W(lpo).next_W(lpo)";
      if(i==1) cout << ".next_W(lpo)";
      if(i==3) cout << ".next_E(lpo)";
      if(i==4) cout << ".next_E(lpo).next_E(lpo)";
      
      if(j==0) cout << ".next_S(lpo).next_S(lpo)";
      if(j==1) cout << ".next_S(lpo)";
      if(j==3) cout << ".next_N(lpo)";
      if(j==4) cout << ".next_N(lpo).next_N(lpo)";
      
      if(k==0) cout << ".next_D(lpo).next_D(lpo)";
      if(k==1) cout << ".next_D(lpo)";
      if(k==3) cout << ".next_T(lpo)";
      if(k==4) cout << ".next_T(lpo).next_T(lpo)";
    }
  }
}

void Print_d_name(dir_3D d) {
  if(d==Wdir) cout << "Wdir";
  if(d==Edir) cout << "Edir";
  if(d==Ndir) cout << "Ndir";
  if(d==Sdir) cout << "Sdir";
  if(d==Tdir) cout << "Tdir";
  if(d==Ddir) cout << "Ddir";
}


void Print_bos_for_restriction(bool idem,
                               dir_3D d,
                               int i_from, int j_from, int k_from,
                               int i_to,   int j_to,   int k_to,
                               int reg_formula[5][5][5]) {
  int value_from, value_to;
  if(i_from >= 0 && i_from < 5 &&
     j_from >= 0 && j_from < 5 &&
     k_from >= 0 && k_from < 5 &&
     
     i_to >= 0 && i_to < 5 &&
     j_to >= 0 && j_to < 5 &&
     k_to >= 0 && k_to < 5) {
    value_from = reg_formula[i_from][j_from][k_from];
    value_to   = reg_formula[i_to][j_to][k_to];
    if(value_to!=0 || value_from!=0) {
      if(idem) {
        if(Drucke_code) {
          if(vom_feinsten_Gitter) {
            cout << "   near_bo_restriction_fine(next);" << endl;
          }
          else {
            cout << " near_bo_restriction_coarse(I";
            Print_ijk_name(i_from,j_from,k_from);
            cout << ");" << endl;
          }
        }
        else {
          cout << " Idem: "; 
          Print_ijk_name(i_from,j_from,k_from);
          cout << endl;
        }
      }
      else {
        if(Drucke_code) {
          cout << "   near_bo_restriction_fine_bo(next,";
            Print_d_name(d);
            cout << ","  << value_from
                 << ","  << value_to
                 << ");" << endl;
        }
        else {
          cout << " bo poi, direction: "; Print(d);
          cout << " from: "; 
          Print_ijk_name(i_from,j_from,k_from);
          cout << " value: " << value_from;
          
          cout << " to: "; 
          Print_ijk_name(i_to,j_to,k_to);
          cout << " value: " << value_to << endl;
        }
      }
    }
  }
}

void Print_restriction_formula() {
  int reg_formula[5][5][5];
  int i,j,k;

  for(i=0;i<5;i++) for(j=0;j<5;j++) for(k=0;k<5;k++) 
    reg_formula[i][j][k] = 0;

  reg_formula[2][2][2] = 2;
  reg_formula[2][3][2] = 1;
  reg_formula[2][1][2] = 1;
  reg_formula[3][2][2] = 1;
  reg_formula[1][2][2] = 1;
  reg_formula[2][2][3] = 1;
  reg_formula[2][2][1] = 1;

  reg_formula[2][3][1] = 1;
  reg_formula[1][3][2] = 1;
  reg_formula[1][2][3] = 1;
  reg_formula[3][2][1] = 1;
  reg_formula[2][1][3] = 1;
  reg_formula[3][1][2] = 1;

  reg_formula[3][1][3] = 1;
  reg_formula[1][3][1] = 1;
  
  for(i=0;i<5;i++) for(j=0;j<5;j++) for(k=0;k<5;k++) {
    if(Drucke_code) {
      cout << " next = I";
      Print_ijk_name(i,j,k);
      cout << ";" << endl;
      cout << " if(gr->Give_type(next) >= interior) {" << endl;
    }
    Print_bos_for_restriction(true,Wdir,
                              i,j,k,
                              i,j,k,
                              reg_formula);
    if(vom_feinsten_Gitter) {
      Print_bos_for_restriction(false,Wdir,
                                i,j,k,
                                i-1,j,k,
                                reg_formula);
      Print_bos_for_restriction(false,Edir,
                                i,j,k,
                                i+1,j,k,
                              reg_formula);
      Print_bos_for_restriction(false,Ndir,
                                i,j,k,
                                i,j+1,k,
                                reg_formula);
      Print_bos_for_restriction(false,Sdir,
                                i,j,k,
                                i,j-1,k,
                                reg_formula);
      Print_bos_for_restriction(false,Tdir,
                                i,j,k,
                                i,j,k+1,
                                reg_formula);
      Print_bos_for_restriction(false,Ddir,
                                i,j,k,
                                i,j,k-1,
                                reg_formula);
    }
    if(Drucke_code) {
      cout << " };" << endl;
    }

    /*
    if(reg_formula[i][j][k]!=0) 
      Print_ijk_name(i,j,k);
    */
  }
}

---