src/expde/formulas/diffop.h

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
//

// ------------------------------------------------------------
// diffop.h
//
// ------------------------------------------------------------

// -----------------------------------------------------------
//    FFFF    OO     RR     M     M   EEEE  L     N     N
//    F      O  O    R R    M M M M   E     L     N N   N
//    FFF    O  O    RR     M  M  M   EEEE  L     N  N  N
//    F      O  O    R R    M     M   E     L     N   N N
//    F       OO     R  R   M     M   EEEE  LLLL  N     N
// -----------------------------------------------------------

// 1.:  lokale Randzellen-Matrizen
// 2.:  Boundary operators


#ifndef DIFFOP_H_
#define DIFFOP_H_
                 

#define XM0 tet->Xm0
#define YM0 tet->Ym0
#define ZM0 tet->Zm0

#define XM1 tet->Xm1
#define YM1 tet->Ym1
#define ZM1 tet->Zm1

#define XM2 tet->Xm2
#define YM2 tet->Ym2
#define ZM2 tet->Zm2

#define XM3 tet->Xm3
#define YM3 tet->Ym3
#define ZM3 tet->Zm3


#define UU0  bocedata->vars[tet->N0()][num_variable]
#define UU1  bocedata->vars[tet->N1()][num_variable]
#define UU2  bocedata->vars[tet->N2()][num_variable]
#define UU3  bocedata->vars[tet->N3()][num_variable]

#define UAA0 U_array[tet->N0()]
#define UAA1 U_array[tet->N1()]
#define UAA2 U_array[tet->N2()]
#define UAA3 U_array[tet->N3()]

#define Coefficient_triangle  (bocedata->vars[tet->N0()][num_coefficient]+bocedata->vars[tet->N1()][num_coefficient]+bocedata->vars[tet->N2()][num_coefficient]+bocedata->vars[tet->N3()][num_coefficient])/4.0

     // 1.: lokale Randzellen-Matrizen

// definiert in boundary.h






   // 2.:  interpolation operators


/*
// interpolation operator for Zellfreiheitsgrad
inline double BoCeData::Interpolation_interior_corners(double *UU) const {


  cout << " weg! " << endl;


  return 1.0;
};

// weight for Zellfreiheitsgrad
inline double BoCeData::Weight_bocellpoint() const {
  int anz, i;

  // arithmetic middle
  anz = 0;
  for(i=0;i<number_points;++i) {
    if(edge_point(i)==corner_poi_typ) {
      ++anz;
    }
  }
  if(developer_version) {
    if(anz==0) cout << "Error in BoCeData::Interpolation_interior_corners! "
                    << endl;
  }
  return 1.0 / anz;
};



//  alless weggggggggggggggggggg::
inline void BoCeData::Nodal_vector_Dirichlet(double *UN, dir_sons i) const {
  int k, anz;
  // Calculate UN
  // a) corner and edge points
  anz = 0;
  for(k=0;k<number_points;++k) { 
    if(edge_point(k) == corner_poi_typ) ++anz;
    if(edge_point(k) == corner_poi_typ && corner(k) == i) {
      UN[k]=1.0;
    }
    else UN[k]=0.0;
  }
  // b) Zellfreiheitsgrad
  // arithmetic middle
  UN[number_points] = 1.0/ anz;
  if(developer_version) {
    if(anz==0) cout << "Error in BoCeData::Nodal_vector_Dirichlet! "
                    << endl;
  }
}

inline void BoCeData::Nodal_vector_Neumann(double *UN, dir_sons i) const {
  int k, anz;
  // Calculate UN
  // a) corner and edge points
  anz=0;
  for(k=0;k<number_points;++k) {
    if(edge_point(k) == corner_poi_typ) ++anz;
    if((corner(k) == i) &&
       (edge_point(k) == corner_poi_typ || edge_point(k) == edge_poi_typ)) {
      UN[k]=1.0;
    }
    else UN[k]=0.0;
  }

  // b) Zellfreiheitsgrad
  // arithmetic middle
  UN[number_points] = 1.0/ anz;
  if(developer_version) {
    if(anz==0) cout << "Error in BoCeData::Nodal_vector_Neumann! "
                    << endl;
  }
}

inline void BoCeData::Nodal_vector_Mixed(double *UN, dir_sons i,
                                         bool* label_on_cell) const {
  int k, anz;
  // Calculate UN
  // a) corner and edge points
  anz=0;
  for(k=0;k<number_points;++k) {
    if(edge_point(k) == corner_poi_typ) ++anz;
    if((corner(k) == i) &&
       (edge_point(k) == corner_poi_typ ||
        (edge_point(k) == edge_poi_typ  && label_on_cell[k] ))) {
      UN[k]=1.0;
    }
    else UN[k]=0.0;
  }

  // b) Zellfreiheitsgrad
  // arithmetic middle
  UN[number_points] = 1.0/ anz;
}
*/

// 2. boundary operators

// L2 Integration
class L2boundary_const {
 public:
  diff_type typ_u() { return  Helm_type; }; 
  diff_type typ_v() { return  Helm_type; };

  L2boundary_const() { }
  static inline int    ops_interior() { return 0; }
  static inline int    ops_diag_interior() { return 0; }
  static inline double stencil(double uN,
                               double uW, double uM, double uO,
                               double uS, double uT, double uD,
                               double uND, double uWN, double uWT,
                               double uED, double uST, double uES,
                               double uEST, double uWND,
                               double h_mesh) {
    return 0.0;
  }
  static inline double diag_stencil(double h_mesh) {
    return 0.0;
  }
  static inline double F_reg_cell(double*const * u_Recell, int num_var,
                                  dir_sons dir_v,
                                  double h_mesh, double* sten) {
    return 0.0;
  }
  static inline double Factor_reg_cell(double h_mesh) {
    return 0.0;
  }
  static inline double diag_F_reg_cell(dir_sons dir_v,
                                       double h_mesh, double* sten) {
    return 0.0;
  }

  // Operator applied to variable with number
  static inline double F_bo_cell(const BoCeData* bocedata,
                                 int num_v, int num_variable) {
    static double sum;
    Tetraeder_storage* tet;
    
    if(bocedata->edge_point(num_v)!=edge_poi_typ) return 0.0;
    else {
      sum=0.0;
      for(tet=bocedata->Give_boundary_tets();tet!=NULL;tet=tet->Next()) {
        if(tet->N0()==num_v) {
          sum = sum + tet->Give_det_surface() *
            (2.0 * UU0 + UU1 + UU2) / 24.0;
        }
        else if(tet->N1()==num_v) {
          sum = sum + tet->Give_det_surface() *
            (2.0 * UU1 + UU0 + UU2) / 24.0;
        }
        else if(tet->N2()==num_v) {
          sum = sum + tet->Give_det_surface() *
            (2.0 * UU2 + UU0 + UU1) / 24.0;
        }
      }  
      return sum;
    }
  }
  static inline double Local_matrix_bo_cell(const BoCeData* bocedata,
                                          double*  U_array, int num_v) {
    static double sum;
    Tetraeder_storage* tet;
    
    if(bocedata->edge_point(num_v)!=edge_poi_typ) return 0.0;
    else {
      sum=0.0;
      for(tet=bocedata->Give_boundary_tets();tet!=NULL;tet=tet->Next()) {
        if(tet->N0()==num_v) {
          sum = sum + tet->Give_det_surface() *
            (2.0 * UAA0 + UAA1 + UAA2) / 24.0;
        }
        else if(tet->N1()==num_v) {
          sum = sum + tet->Give_det_surface() *
            (2.0 * UAA1 + UAA0 + UAA2) / 24.0;
        }
        else if(tet->N2()==num_v) {
          sum = sum + tet->Give_det_surface() *
            (2.0 * UAA2 + UAA0 + UAA1) / 24.0;
        }
      }  
      return sum;
    }
  }
  static inline double diag_F_bo_cell(const BoCeData* bocedata,
                                      int num_v) {
    static double sum;
    Tetraeder_storage* tet;
    
    if(bocedata->edge_point(num_v)!=edge_poi_typ) return 0.0;
    else {
      sum=0.0;
      for(tet=bocedata->Give_boundary_tets();tet!=NULL;tet=tet->Next()) {
        if(tet->N0()==num_v) {
          sum = sum + tet->Give_det_surface() / 12.0;
        }
        else if(tet->N1()==num_v) {
          sum = sum + tet->Give_det_surface() / 12.0;
        }
        else if(tet->N2()==num_v) {
          sum = sum + tet->Give_det_surface() / 12.0;
        }
      }  
      return sum;
    }
  }
};


// L2 Integration variable
class L2boundary_variable {
 public:
  diff_type typ_u() { return  Helm_type; }; 
  diff_type typ_v() { return  Helm_type; };

  L2boundary_variable() { }
  static inline int    ops_interior() { return 0; }
  static inline int    ops_diag_interior() { return 0; }
  static inline double stencil(double uN,
                               double uW, double uM, double uE,
                               double uS, double uT, double uD,
                               double uND, double uWN, double uWT,
                               double uED, double uST, double uES,
                               double uEST, double uWND,
                               double aN,
                               double aW, double aM, double aE,
                               double aS, double aT, double aD,
                               double aST, double aET, double aES,
                               double aND, double aWD, double aWN,
                               double aEST, double aWND,
                               double h_mesh) {
    return 0.0;
  }
  static inline double diag_stencil(double h_mesh,
                                    double aN,
                                    double aW, double aM, double aE,
                                    double aS, double aT, double aD,
                                    double aST, double aET, double aES,
                                    double aND, double aWD, double aWN,
                                    double aEST, double aWND) {
    return 0.0;
  }
  static inline double F_reg_cell(double*const * u_Recell, int num_var,
                                  dir_sons dir_v,
                                  double h_mesh, double* sten) {
    return 0.0;
  }
  static inline double Factor_reg_cell(double h_mesh) {
    return 0.0;
  }
  static inline double diag_F_reg_cell(dir_sons dir_v,
                                       double h_mesh, double* sten) {
    return 0.0;
  }

  // Operator applied to variable with number
  static inline double F_bo_cell(const BoCeData* bocedata,
                                 int num_v, int num_variable,
                                 int num_coefficient) {
    static double sum;
    Tetraeder_storage* tet;
    
    if(bocedata->edge_point(num_v)!=edge_poi_typ) return 0.0;
    else {
      sum=0.0;
      for(tet=bocedata->Give_boundary_tets();tet!=NULL;tet=tet->Next()) {
        if(tet->N0()==num_v) {
          sum = sum + tet->Give_det_surface() *
            (2.0 * UU0 + UU1 + UU2) / 24.0 * Coefficient_triangle;
        }
        else if(tet->N1()==num_v) {
          sum = sum + tet->Give_det_surface() *
            (2.0 * UU1 + UU0 + UU2) / 24.0 * Coefficient_triangle;
        }
        else if(tet->N2()==num_v) {
          sum = sum + tet->Give_det_surface() *
            (2.0 * UU2 + UU0 + UU1) / 24.0 * Coefficient_triangle;
        }
      }  
      return sum;
    }
  }
  static inline double Local_matrix_bo_cell(const BoCeData* bocedata,
                                            double*  U_array, int num_v,
                                            int num_coefficient) {
    static double sum;
    Tetraeder_storage* tet;
    
    if(bocedata->edge_point(num_v)!=edge_poi_typ) return 0.0;
    else {
      sum=0.0;
      for(tet=bocedata->Give_boundary_tets();tet!=NULL;tet=tet->Next()) {
        if(tet->N0()==num_v) {
          sum = sum + tet->Give_det_surface() *
            (2.0 * UAA0 + UAA1 + UAA2) / 24.0 * Coefficient_triangle;
        }
        else if(tet->N1()==num_v) {
          sum = sum + tet->Give_det_surface() *
            (2.0 * UAA1 + UAA0 + UAA2) / 24.0 * Coefficient_triangle;
        }
        else if(tet->N2()==num_v) {
          sum = sum + tet->Give_det_surface() *
            (2.0 * UAA2 + UAA0 + UAA1) / 24.0 * Coefficient_triangle;
        }
      }  
      return sum;
    }
  }
  static inline double diag_F_bo_cell(const BoCeData* bocedata,
                                      int num_v, int num_coefficient) {
    static double sum;
    Tetraeder_storage* tet;
    
    if(bocedata->edge_point(num_v)!=edge_poi_typ) return 0.0;
    else {
      sum=0.0;
      for(tet=bocedata->Give_boundary_tets();tet!=NULL;tet=tet->Next()) {
        if(tet->N0()==num_v) {
          sum = sum + tet->Give_det_surface() / 12.0 * Coefficient_triangle;
        }
        else if(tet->N1()==num_v) {
          sum = sum + tet->Give_det_surface() / 12.0 * Coefficient_triangle;
        }
        else if(tet->N2()==num_v) {
          sum = sum + tet->Give_det_surface() / 12.0 * Coefficient_triangle;
        }
      }  
      return sum;
    }
  }
};

#endif
        

---