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

// ------------------------------------------------------------
// variable.cc
//
// ------------------------------------------------------------

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



// Include level 0:
#include "../parser.h"

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

// Include level 2:
#include "../basic/basic.h"

// Include level 3:
#include "../domain/domain.h"

// Include level 4:
#include "../formulas/boundy.h"
#include "../formulas/loc_sten.h"

// Include level 5:
#include "../grid/gpar.h"
#include "../grid/parallel.h"
#include "../grid/mgcoeff.h"
#include "../grid/sto_man.h"
#include "../grid/gridbase.h"
#include "../grid/grid.h"
#include "../grid/input.h"


// Include level 5.1
#include "../evpar/evpar.h"

// Include level 6:
#include "../extemp/variable.h"  
#include "../extemp/eleme.h"
     
// Include level 7:
#include "index.h"


               // 1. Index_set
               // 2. ExpdeIndex
               // 3. from variable.h

// 1. Index_set
Index_set::Index_set(Grid* gr) {
  int max_size;

  // -1.
  hashtable_initialized = false;
  nicht_implementiert   = false;

  // 0. 
  grid = gr;

  // 2. construct elements
  used_P_interior = 0;
  used_P_nearb    = 0;
  used_P_cellpoi  = 0;
  used_P_Bo2p     = 0;

  length_P_interior =  minimal_length_P_interior;
  length_P_nearb    =  minimal_length_P_nearb;
  length_P_cellpoi  =  minimal_length_P_cellpoi;
  length_P_Bo2p     =  minimal_length_P_Bo2p;

  start_P_interior  = new P_interior*[length_P_interior];
  start_P_nearb     = new P_nearb*[length_P_nearb];
  start_P_cellpoi   = new P_cellpoi*[length_P_cellpoi];
  start_P_Bo2p      = new P_Bo2p*[length_P_Bo2p];

  size     = used_P_interior + used_P_nearb + used_P_cellpoi + used_P_Bo2p;
  max_size = 
    length_P_interior + length_P_nearb + 
    length_P_cellpoi  + length_P_Bo2p;

  // sorted numbers
  sorted_number = new int[max_size];
  sorted_number_back = new int[max_size];
  int ii;
  for(ii=0;ii<max_size;++ii) {
    sorted_number[ii] = ii;
    sorted_number_back[ii] = ii;
  }
}



Index_set::Index_set(Subgrid subgrid, Grid* gr) {
  int color;
  int lev;
  int run_bo;
  int number;
 
  // -1.
  hashtable_initialized = false;
  nicht_implementiert   = false;

  sorted_number = NULL;

  // 0. 
  P_interior *iter_i;
  P_nearb    *iter_n;
  P_Bo2p     *iter_b;
  P_cellpoi  *iter_cf;


  grid = gr;
  lev = grid->Max_level();

  // 1. count
  used_P_interior = 0;
  if(subgrid.run_interior()) {
    // innere Punkte
    for(color=0;color<8;++color)
      for(iter_i=grid->Start_P_interior(lev,color);
          iter_i!=NULL;iter_i=iter_i->Next()) {
        ++used_P_interior;
      }
  }
  used_P_nearb = 0;
  used_P_cellpoi = 0;
  if(subgrid.run_nearb()) {
    // Randzellfreiheitsgrade
    for(iter_cf=grid->Start_P_cellpoi(lev);
        iter_cf!=NULL;iter_cf=iter_cf->Next()) {
      ++used_P_cellpoi;
    }
    // randnahe Punkte
    for(color=0;color<8;++color)
      for(iter_n=grid->Start_P_nearb(lev,color);
            iter_n!=NULL;iter_n=iter_n->Next()) {
        ++used_P_nearb;
      }
  }
  used_P_Bo2p = 0;
  run_bo = subgrid.run_boundary();
  if(run_bo==true) {
    // Randpunkte
    for(color=0;color<8;++color) {
      for(iter_b=grid->Start_P_Bo2p(lev,color);
          iter_b!=NULL;iter_b=iter_b->Next()) {
        ++used_P_Bo2p;
      }
    }
  }
  else if(run_bo<false) {
    // Randpunkte
    for(color=0;color<8;++color)
      for(iter_b=grid->Start_P_Bo2p(lev,color);
          iter_b!=NULL;iter_b=iter_b->Next()) {
        if(iter_b->Give_Label(-run_bo,grid)) {
          ++used_P_Bo2p;
        }
      }
  }
  
  // 2. construct elements
  length_P_interior = used_P_interior;
  length_P_nearb    = used_P_nearb;
  length_P_cellpoi  = used_P_cellpoi;
  length_P_Bo2p     = used_P_Bo2p;

  if(length_P_interior < minimal_length_P_interior) 
    length_P_interior =  minimal_length_P_interior;
  if(length_P_nearb    < minimal_length_P_nearb)    
    length_P_nearb    =  minimal_length_P_nearb;
  if(length_P_cellpoi  < minimal_length_P_cellpoi)
    length_P_cellpoi  =  minimal_length_P_cellpoi;
  if(length_P_Bo2p     < minimal_length_P_Bo2p)     
    length_P_Bo2p     =  minimal_length_P_Bo2p;

  start_P_interior = new P_interior*[length_P_interior];
  start_P_nearb    = new P_nearb*[length_P_nearb];
  start_P_cellpoi  = new P_cellpoi*[length_P_cellpoi];
  start_P_Bo2p     = new P_Bo2p*[length_P_Bo2p];


  size = used_P_interior + used_P_nearb + used_P_cellpoi + used_P_Bo2p;

  // 3. set pointers
  if(subgrid.run_interior()) {
    number = 0;
    // innere Punkte
    for(color=0;color<8;++color)
      for(iter_i=grid->Start_P_interior(lev,color);
          iter_i!=NULL;iter_i=iter_i->Next()) {
        start_P_interior[number] = iter_i;
        ++number;
      }
  }

  if(subgrid.run_nearb()) {
    number=0;
    // Randzellfreiheitsgrade
    for(iter_cf=grid->Start_P_cellpoi(lev);
        iter_cf!=NULL;iter_cf=iter_cf->Next()) {
        start_P_cellpoi[number] = iter_cf;
        ++number;
    }
    number=0;
    // randnahe Punkte
    for(color=0;color<8;++color)
      for(iter_n=grid->Start_P_nearb(lev,color);
            iter_n!=NULL;iter_n=iter_n->Next()) {
        start_P_nearb[number] = iter_n;
        ++number;
      }
  }

  run_bo = subgrid.run_boundary();
  if(run_bo==true) {
    number = 0;
    // Randpunkte
    for(color=0;color<8;++color) {
      for(iter_b=grid->Start_P_Bo2p(lev,color);
          iter_b!=NULL;iter_b=iter_b->Next()) {
        start_P_Bo2p[number] = iter_b;
        ++number;
      }
    }
  }
  else if(run_bo<false) {
    number = 0;
    // Randpunkte
    for(color=0;color<8;++color)
      for(iter_b=grid->Start_P_Bo2p(lev,color);
          iter_b!=NULL;iter_b=iter_b->Next()) {
        if(iter_b->Give_Label(-run_bo,grid)) {
        start_P_Bo2p[number] = iter_b;
        ++number;
        }
      }
  }

  // sorted numbers
  sorted_number = new int[size];
  sorted_number_back = new int[size];
  int ii;
  for(ii=0;ii<size;++ii) {
    sorted_number[ii] = ii;
    sorted_number_back[ii] = ii;
  }
};

// for sort_planes
int Index_set::integer_m(int ii,
                         double coordz_min, 
                         int num_points_z,
                         int zwei_po_lev) {
  typ_of_general_points typ;
  Index3D I;
  dir_3D d;

  typ = type(ii);
  if(typ==type_P_Bo2p) {
    // Bo2p
    d = pointer_P_Bo2p(ii)->d();
    if(d==Ddir) return 0;
    if(d==Tdir) return num_points_z-1;

    I = pointer_P_Bo2p(ii)->Ind();
    return (int)((I.coordinate_z()-coordz_min) * zwei_po_lev) + 1;
  }
  else {
    if(typ==type_P_interior) {
      // interior
      I = pointer_P_interior(ii)->Ind();
    }
    if(typ==type_P_nearb) {
      // nearb
      I = pointer_P_nearb(ii)->Ind();
    }
    return (int)((I.coordinate_z()-coordz_min) * zwei_po_lev) + 1;
  }
};


void Index_set::Sort_by(sorting_types sort_typ) {
  int m;
  int zwei_po_lev;

  Index3D I;
  double coordz_max, coordz_min;
  ExpdeIndex i;
  int ii;
  int num_points_z;
  int *anz_in_plane;

  Index1D Iz_min;
  Index1D Iz_max;

  typ_of_general_points typ;
  
  if(sort_typ != z_plane_sort) {
    cout << " Dieser Sortierungstyp ist in Index_set::Index_set " 
         << " nicht implementiert! " << endl;
  }

  zwei_po_lev = Zweipotenz(grid->Max_level());

  // 1. count and look for small and large index3D
  coordz_max = 0.0;
  coordz_min = 1.0;
 
  for(ii=0;ii<size;++ii) {
    typ = type(ii);
    if(typ==type_P_interior) {
      // interior
      I = pointer_P_interior(ii)->Ind();
    }
    if(typ==type_P_nearb) {
      // nearb
      I = pointer_P_nearb(ii)->Ind();
    }
    if(typ==type_P_cellpoi) {
      cout << " Noch nicht implementiert in diesem Index Konstruktor! "
           << endl;
    }
    if(typ==type_P_Bo2p) {
      // Bo2p
      I = pointer_P_Bo2p(ii)->Ind();
    }
    // maximal and minimal point in z-direction
    if(coordz_max < I.coordinate_z()) {
      Iz_max     = I.I_z();
      coordz_max = I.coordinate_z();
    }
    if(coordz_min > I.coordinate_z()) {
      Iz_min     = I.I_z();
      coordz_min = I.coordinate_z();
    }
  }

  // 3. count number of points in z direction
  num_points_z = (int)((coordz_max-coordz_min) *
                       zwei_po_lev) + 3;
  anz_in_plane = new int[num_points_z];
  for(m=0; m<num_points_z; ++m) {
    anz_in_plane[m] = 0;
  }

  // 4. count anz_in_plane[m]
  for(ii=0;ii<size;++ii) {
    m  = integer_m(ii,coordz_min,num_points_z,zwei_po_lev);
    anz_in_plane[m] = anz_in_plane[m] + 1;
  }

  for(m=1; m<num_points_z; ++m) {
    anz_in_plane[m] = anz_in_plane[m] + anz_in_plane[m-1];
  }

  for(m=num_points_z-1; m>0; --m) {
    anz_in_plane[m] = anz_in_plane[m-1];
  }
  anz_in_plane[0] = 0;

  /*
  for(m=0; m<num_points_z; ++m) {
    cout << " anz_in_plane[" << m << "] = " 
         << anz_in_plane[m] << endl;
  }
  */

  // 5. sorted numbers 
  for(ii=0;ii<size;++ii) {
    m  = integer_m(ii,coordz_min,num_points_z,zwei_po_lev);
    sorted_number[anz_in_plane[m]] = ii;
    sorted_number_back[ii] = anz_in_plane[m];
    anz_in_plane[m]++;
  }

  /*
  for(ii=0;ii<size;++ii) {
    cout << " sorted_number[" << ii << "] = " 
         << sorted_number[ii] << endl;
  }
  */

  /*
  int sum, sumback, sumtest;
  sum     = 0;
  sumback = 0;
  sumtest = 0;

  for(ii=0;ii<size;++ii) {
    sum     = sum     + sorted_number[ii];
    sumtest = sumtest + ii;
    sumback = sumback + sorted_number_back[ii];
  }

  cout << " sum: "     << sum 
       << " sumtest: " << sumtest 
       << " sumback: " << sumback << endl;
  */
};

void Index_set::initialize_hashtable() {
  int i;

  if(nicht_implementiert) 
    cout << " error in Index_set::initialize_hashtable() nicht implementiert"
         << endl;

  if(hashtable_initialized==false) {
    hashtable_initialized = true;
    // hashtable initialisieren
    hashtable_indices_leng  = 
      Find_next_prime(Factor_hashtable_indices_lenght*size);
    hashtable_indices_start = 
      new Point_hashtable_indices*[hashtable_indices_leng];
    for(i=0;i<hashtable_indices_leng;++i)
      hashtable_indices_start[i] = NULL;
    hashtable_indices_occ=0;
    // Werte setzen
    for(i=0;i<used_P_interior;++i) {
      Add_point(start_P_interior[i]->Ind(),6,i);
    }
    for(i=0;i<used_P_nearb;++i) {
      Add_point(start_P_nearb[i]->Ind(),6,i+used_P_interior);
    }
    for(i=0;i<used_P_cellpoi;++i) {
      Add_point(start_P_cellpoi[i]->Ind(),7,
                i+used_P_interior+used_P_nearb);
    }
    for(i=0;i<used_P_Bo2p;++i) {

      Add_point(start_P_Bo2p[i]->Ind(),(int)start_P_Bo2p[i]->d(),
                i+used_P_interior+used_P_nearb+used_P_cellpoi);
    }
  }
}


void Index_set::Add_point(Index3D I, int Gedir, int value) {
  int Nx, Ny, Nz;
  Nx = I.I_x().get();
  Ny = I.I_y().get();
  Nz = I.I_z().get();

  Point_hashtable_indices* point;

  point = hashtable_indices_start
    [hashtable_indices_function(Nx,Ny,Nz,Gedir,hashtable_indices_leng)];
  if(point==NULL) {
    hashtable_indices_occ++;
    point = new Point_hashtable_indices(Nx,Ny,Nz,Gedir,value);
    hashtable_indices_start[hashtable_indices_function(Nx,Ny,Nz,Gedir,
                                                   hashtable_indices_leng)] = 
      point;
  }
  else {
    while(point->next!=NULL && 
          (point->nx!=Nx || point->ny!=Ny || point->nz!=Nz
           || point->gedir!=Gedir))
      point = point->next;
    if(point->nx!=Nx || point->ny!=Ny || point->nz!=Nz
       || point->gedir!=Gedir) {
      hashtable_indices_occ++;
      point->next = new Point_hashtable_indices(Nx,Ny,Nz,Gedir,value);
    }
  }
}


int Index_set::Give_value(Index3D I, int Gedir) {
  int Nx, Ny, Nz;
  Nx = I.I_x().get();
  Ny = I.I_y().get();
  Nz = I.I_z().get();

  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->wert;
    point = point->next;
  }
  cout << " Error in Index_set::Give_value(Index3D I, int Gedir) " << endl;
  return 0;      
}

bool Index_set::Exists_index(Index3D I, int Gedir) {
  int Nx, Ny, Nz;
  Nx = I.I_x().get();
  Ny = I.I_y().get();
  Nz = I.I_z().get();

  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 true;
    point = point->next;
  }
  return false;      
}


// 2. ExpdeIndex

ExpdeIndex::ExpdeIndex() {
  set = NULL;
  i = 0;
  counter = 0;
}

ExpdeIndex::ExpdeIndex(int co, Index_set* s) {
  set     = s;
  counter = co;
  i       = set->Sorted_number(counter);
}


void ExpdeIndex::Test_correct() {
  if(set == NULL) cout << " Error: set == NULL in ExpdeIndex! " << endl;
  else {
    if(i < 0) cout << " Error: i < 0 in ExpdeIndex! " << endl;
    if(i >= set->Size()) cout << " Error: i >= Size in ExpdeIndex! " << endl;
  }
}

int ExpdeIndex::integer(Index_set& M) {
  int gedir;
  Index3D I;
  typ_of_general_points typ;
  
  // 0. short cut?
  if(&M == set) {
    //    return i;
    return counter;
  }
  // 1. calc type 
  typ = set->type(i);

  // calc Index 3D and gedir;
  if(typ==type_P_interior) {
    I = set->pointer_P_interior(i)->Ind();
    gedir = 6;
  }
  else if(typ==type_P_nearb){
    I = set->pointer_P_nearb(i)->Ind();
    gedir = 6;
  } 
  else if(typ==type_P_Bo2p){
    I = set->pointer_P_Bo2p(i)->Ind();
    gedir = set->pointer_P_Bo2p(i)->d();
  }
  else if(typ==type_P_cellpoi){
    I = set->pointer_P_cellpoi(i)->Ind();
    gedir = 7;
  }

  return M.Sorted_number_back(M.Give_value(I,gedir));
  //  return M.Give_value(I,gedir);
}


// 3. from variable.h

assignable_vector Variable::operator[] (ExpdeIndex& i) {
  return assignable_vector(i.index_set()->varM(i.i_integer()),
                           number_variable);
}

---