src/expde/indices/index.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
//

// ------------------------------------------------------------
//
//     index.h 
//
// ------------------------------------------------------------

#ifndef INDEX_H_
#define INDEX_H_
       
//  Index_set, ExpdeIndex

class Index_set;
class ExpdeIndex;

#include "grid/SSten.h"

// sorting types
enum sorting_types { z_plane_sort };

// a) hash table data structure

class Point_hashtable_indices {
  friend class Index_set;
 private:
  int nx, ny, nz, gedir;   // inter code

  int wert;

  Point_hashtable_indices* next;
 public:
  Point_hashtable_indices(int Nx, int Ny, int Nz, int Gedir, 
                          int value) 
    : nx(Nx), ny(Ny), nz(Nz), gedir(Gedir), wert(value) {
    next = NULL;
  };
  int Value() { return wert; }
  ~Point_hashtable_indices() {
    delete next;
  };
};

inline int hashtable_indices_function(int x, int y, int z, int d, int length) {
  return (d + 11*x + 101*y + 10007*z) % length;
}




// b) ...


#define minimal_length_P_interior  27
#define minimal_length_P_nearb     27
#define minimal_length_P_cellpoi   30
#define minimal_length_P_Bo2p      20


class assignable_vector {
 private:
  double *var;
  int num;
 public:
  assignable_vector(double* v, int nu) : var(v), num(nu) {};
  void operator=(double w) { var[num]=w; };
};

class Index_set_Expr {
 public:
  Index_set* set;
  ExpdeIndex* i;
  Index_set_Expr(Index_set* Kp, ExpdeIndex* ii) : set(Kp), i(ii) {};
};

class Index_set {
friend class Variable;
 public:
  Index_set(Subgrid subgrid, Grid* gr);
  Index_set(Grid* gr);

  void Sort_by(sorting_types sort_typ);

  int Size() { return size; }
  ExpdeIndex first();
  ExpdeIndex last();
  Grid* Give_grid() { return grid; }

  void operator=(Index_set_Expr index_set_expr);

  typ_of_general_points type(int n);

  P_interior *pointer_P_interior(int n);
  P_nearb    *pointer_P_nearb(int n);
  P_cellpoi  *pointer_P_cellpoi(int n);
  P_Bo2p     *pointer_P_Bo2p(int n);

  // give number i of element a_i in this set
  int  Give_value(Index3D I, int Gedir);
  int  Sorted_number(int i) { return sorted_number[i]; }
  int  Sorted_number_back(int i) { return sorted_number_back[i]; }
 private:
  // give pointer to doubles in grid
  double *varM(int l);

  // grid
  Grid* grid;

  // vector data structure
  int size;
  // numbering:
  // 1. P_interior
  // 2. P_nearb
  // 3. P_cellpoi
  // 4. P_Bo2p

  // for sort_planes
  int *sorted_number;
  int *sorted_number_back;
  int integer_m(int ii, double coordz_min,
                int num_points_z,
                int zwei_po_lev);


  // starting points of elements
  P_interior **start_P_interior;
  P_nearb    **start_P_nearb;
  P_cellpoi  **start_P_cellpoi;
  P_Bo2p     **start_P_Bo2p;

  // length
  int length_P_interior;
  int length_P_nearb;
  int length_P_cellpoi;
  int length_P_Bo2p;

  // used
  int used_P_interior;
  int used_P_nearb;
  int used_P_cellpoi;
  int used_P_Bo2p;

  // set new array's
  void Start_setting_arrays();

  void Add_regular_condi(Index3D I, Index_set* set);
  void Add_cellpoi_condi(Index3D I, Index_set* set);
  void Add_Bo2p_condi(Index3D I, dir_3D d, Index_set* set);

  // fuer hashtable
  // Punkt hinzufuegen usw.
  void Add_point(Index3D I, int Gedir, int value);
  bool Exists_index(Index3D I, int Gedir);

  // hashtable data structure
  bool hashtable_initialized;
  bool nicht_implementiert;
  void initialize_hashtable();

  Point_hashtable_indices* point;              // fuer Iterationen
  int hashtable_indices_leng;                  // Laenge des hashtable
  int hashtable_indices_occ;                   // Besetzungszahl des hashtable
  Point_hashtable_indices** hashtable_indices_start;   // Zeiger auf Tabelle

  // gehe zu hashtable
  Point_hashtable_indices* hashtable_indices_point(int Nx, int Ny, int Nz,
                                                   int Gedir) const;

  // Abkuerzung fuer  gedir:
  // Bo2p:            gedir = 0,...,5  dir_3D
  // interior, nearb: gedir = 6
  // cell_poi:        gedir = 7


};


class ExpdeIndex {
 private:
  int counter;
  int i;
  Index_set *set;
  void Test_correct();
 public:
  ExpdeIndex();
  ExpdeIndex(int ii, Index_set* s);

  int i_integer() { 
    if(developer_version) Test_correct();
    return i; 
  };
  int integer() { 
    if(developer_version) Test_correct();
    //return i; 
    return counter; 
  };
  int integer(Index_set& M);
  Index_set* index_set() {
    if(developer_version) Test_correct();
    return set;
  }
  void operator++();
  void operator--();
  Index_set_Expr neighbours_in(Index_set& K);
  bool operator<=(const ExpdeIndex& IR);
  bool operator>=(const ExpdeIndex& IR);
  Grid* Give_grid() { return set->Give_grid(); }
};



// some inline functions

inline Point_hashtable_indices* 
Index_set::hashtable_indices_point(int Nx, int Ny, int Nz,
                                   int Gedir) const {
  Point_hashtable_indices* point;
  point = hashtable_indices_start[hashtable_indices_function(Nx,Ny,Nz,Gedir,
                                 hashtable_indices_leng)];
  while(point!=NULL) {
    if(point->nx==Nx && point->ny==Ny && point->nz==Nz &&
       point->gedir==Gedir) return point;
    point = point->next;
  }
  return NULL;
}



inline ExpdeIndex Index_set::first() {
  return ExpdeIndex(0,this);
};

inline ExpdeIndex Index_set::last() {
  return ExpdeIndex(size-1,this);
};


  // give pointer to doubles in grid
inline double *Index_set::varM(int n) {
  if(n < used_P_interior) 
    return start_P_interior[n]->varM(grid,grid->Max_level());
  if(n < used_P_nearb+used_P_interior)
    return start_P_nearb[n-used_P_interior]->varM(grid,grid->Max_level());
  if(n < used_P_cellpoi+used_P_nearb+used_P_interior) 
    return start_P_cellpoi[n-used_P_nearb-used_P_interior]->varM(grid);
  if(developer_version) {
    if(n>=size) 
      cout << " error in Index_set::varM(int n); " 
           << endl;
  }
  return start_P_Bo2p[n-used_P_nearb-used_P_interior
                     -used_P_cellpoi]->varM(grid);
}

 // starting points of elements
inline P_interior *Index_set::pointer_P_interior(int n) {
  if(developer_version) {
    if(n<0 || n>=used_P_interior) 
      cout << " error in Index_set::pointer_P_interior " << endl;
  }
  return start_P_interior[n];
}

inline P_nearb *Index_set::pointer_P_nearb(int n) {
  if(developer_version) {
    if(n<used_P_interior || n>=used_P_nearb+used_P_interior) 
      cout << " error in Index_set::pointer_P_nearb " << endl;
  }
  return start_P_nearb[n-used_P_interior];
}

inline P_cellpoi *Index_set::pointer_P_cellpoi(int n) {
  if(developer_version) {
    if(n<used_P_interior+used_P_nearb || 
       n>=used_P_nearb+used_P_interior+used_P_cellpoi) 
      cout << " error in Index_set::pointer_P_cellpoi " << endl;
  }
  return start_P_cellpoi[n-used_P_interior-used_P_nearb];
}

inline P_Bo2p *Index_set::pointer_P_Bo2p(int n) {
  if(developer_version) {
    if(n<used_P_interior+used_P_nearb+used_P_cellpoi || 
       n>=size) 
      cout << " error in Index_set::pointer_P_Bo2p " 
           << n 
           << " u: " << used_P_interior+used_P_nearb+used_P_cellpoi
           << " o: " << size << endl;
  }
  return start_P_Bo2p[n-used_P_interior-used_P_nearb-used_P_cellpoi];
}



inline typ_of_general_points Index_set::type(int n) {
  if(n < used_P_interior)
    return type_P_interior;
  if(n < used_P_nearb+used_P_interior)
    return type_P_nearb;
  if(n < used_P_cellpoi+ used_P_nearb+used_P_interior)
    return type_P_cellpoi;
  if(developer_version) {
    if(n>=size) 
      cout << " error in Index_set::type(int n); " 
           << endl;
  }
  return type_P_Bo2p;
}

inline Index_set_Expr ExpdeIndex::neighbours_in(Index_set& K) {
  return Index_set_Expr(&K,this);
}


inline void ExpdeIndex::operator++() { 
  if(developer_version) Test_correct();

  ++counter; 
  i = set->Sorted_number(counter);
};

inline void ExpdeIndex::operator--() { 
  if(developer_version) Test_correct();

  --counter; 
  i = set->Sorted_number(counter);
};

inline bool ExpdeIndex::operator<=(const ExpdeIndex& IR) {
  // cout << " Vergl: " << i << ", and " << IR.i << endl;
  return counter <= IR.counter;
}
inline bool ExpdeIndex::operator>=(const ExpdeIndex& IR) {
  // cout << " Vergl: " << i << ", and " << IR.i << endl;
  return counter >= IR.counter;
}




template<class A>
double DExpr<A>::operator[] (ExpdeIndex& i) {
 
  P_interior *iter_i;
  P_nearb    *iter_n;
  P_Bo2p     *iter_b;
  P_cellpoi  *iter_cf;

  typ_of_general_points typ;

  int lev;
  Grid *grid;
  double h_mesh;
  int size_stencil;

  // Zwischenspeicher
  Nearb_Ablage* nearb_ablage;

  grid = i.Give_grid();
  lev  = grid->Max_level();
  h_mesh = grid->H_mesh() / Zweipotenz(lev);

  typ = i.index_set()->type(i.i_integer());
  if(typ==type_P_interior) {
    iter_i = i.index_set()->pointer_P_interior(i.i_integer());
    
    if(developer_version) {
      if(iter_i==NULL) cout << " Error: iter_i " << endl;
    }
    size_stencil=a_.Sice_stencil();
    if(size_stencil<=15)
      return a_.Give_interior(iter_i,grid,h_mesh,lev,
                              iter_i->varN(grid,lev),
                              iter_i->varW(grid,lev),
                              iter_i->varM(grid,lev),
                              iter_i->varE(grid,lev),
                              iter_i->varS(grid,lev),
                              iter_i->varT(grid,lev),
                              iter_i->varD(grid,lev),
                              iter_i->varND(grid,lev),
                              iter_i->varWN(grid,lev),
                              iter_i->varWT(grid,lev),
                              iter_i->varED(grid,lev),
                              iter_i->varST(grid,lev),
                              iter_i->varES(grid,lev),
                              NULL, NULL,
                              iter_i->varEST(grid,lev),
                              iter_i->varWND(grid,lev),
                              NULL, NULL, NULL, NULL,
                              NULL, NULL, NULL, NULL);
    if(size_stencil==17)
      return a_.Give_interior(iter_i,grid,h_mesh,lev,
                              iter_i->varN(grid,lev),
                              iter_i->varW(grid,lev),
                              iter_i->varM(grid,lev),
                              iter_i->varE(grid,lev),
                              iter_i->varS(grid,lev),
                              iter_i->varT(grid,lev),
                              iter_i->varD(grid,lev),
                              iter_i->varND(grid,lev),
                              iter_i->varWN(grid,lev),
                              iter_i->varWT(grid,lev),
                              iter_i->varED(grid,lev),
                              iter_i->varST(grid,lev),
                              iter_i->varES(grid,lev),
                              iter_i->varWD(grid,lev),
                              iter_i->varET(grid,lev),
                              iter_i->varEST(grid,lev),
                              iter_i->varWND(grid,lev),
                              NULL, NULL, NULL, NULL,
                              NULL, NULL, NULL, NULL);
    if(size_stencil==25)
      return a_.Give_interior(iter_i,grid,h_mesh,lev,
                              iter_i->varN(grid,lev),
                              iter_i->varW(grid,lev),
                              iter_i->varM(grid,lev),
                              iter_i->varE(grid,lev),
                              iter_i->varS(grid,lev),
                              iter_i->varT(grid,lev),
                              iter_i->varD(grid,lev),
                              iter_i->varND(grid,lev),
                              iter_i->varWN(grid,lev),
                              iter_i->varWT(grid,lev),
                              iter_i->varED(grid,lev),
                              iter_i->varST(grid,lev),
                              iter_i->varES(grid,lev),
                              iter_i->varWD(grid,lev),
                              iter_i->varET(grid,lev),
                              iter_i->varEST(grid,lev),
                              iter_i->varWND(grid,lev),
                              iter_i->cell_varWSD(grid,lev),
                              iter_i->cell_varWST(grid,lev),
                              iter_i->cell_varWND(grid,lev),
                              iter_i->cell_varWNT(grid,lev),
                              iter_i->cell_varESD(grid,lev),
                              iter_i->cell_varEST(grid,lev),
                              iter_i->cell_varEND(grid,lev),
                              iter_i->cell_varENT(grid,lev));
  }
  
  if(typ==type_P_cellpoi) {
    iter_cf = i.index_set()->pointer_P_cellpoi(i.i_integer());
    
    if(developer_version) {
      if(iter_cf==NULL) cout << " Error: iter_cf " << endl;
    }
    return
      a_.Give_cellpoi(iter_cf,grid,grid->Give_Bo_cell(iter_cf->Ind()));
  }

  if(typ==type_P_nearb) {
    iter_n = i.index_set()->pointer_P_nearb(i.i_integer());

    if(developer_version) {
      if(iter_n==NULL) cout << " Error: iter_n " << endl;
    }
    nearb_ablage = grid->Give_Nearb_Ablage();
    size_stencil=a_.Sice_stencil();
    if(size_stencil>=17) {
      nearb_ablage->Initialize_17(grid,iter_n,lev);
    }
    else {
      nearb_ablage->Initialize(grid,iter_n,lev);
    }
    return
      a_.Give_nearb(iter_n,grid,h_mesh,lev,nearb_ablage);
  }

  if(typ==type_P_Bo2p) {
    iter_b = i.index_set()->pointer_P_Bo2p(i.i_integer());

    if(developer_version) {
      if(iter_b==NULL) cout << " Error: iter_b " << endl;
    }
    return
      a_.Give_Bo2p(iter_b,grid,lev);
  }
  cout << " Error in DExpr<A>::operator[] " << endl;


  cout << " Hello, we are returning zero!!" << endl;
  return 0.0;
};


#endif
         

---