src/expde/grid/mgcoeff.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
//

// ------------------------------------------------------------
//
// mgcoeff.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 "../domain/d_exam.h"

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

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

//  1. Calc_multigrid_points_part2()
//  2. von alter version; muss mal erneuert werden

//  1. Calc_multigrid_points_part2()


void Grid_base::Calc_multigrid_points_part2() {
  int level;
  Index3D I, Iop, Iproc;
  Pointtype typI;
  int i;

  //Test_just_this();


  D3vector v;
  D3vector w;

  if(I_am_active()) {
    // first part: finest points
    // -----------
    level = max_level;
    // interior points
    iterate_hash1 {
      I = point1->Give_Index();
      typI = Give_type(I);
      if(typI>=interior && I.Tiefe()==level) 
        Add_interpolation_point(I,level);
    }
    
    // Zellfreiheitgrade
    iterate_hash2 {
      if(bocell->Exists_bocellpoint()) {
        I = bocell->Give_Index();
        for(i=0;i<8;++i)
          Add_interpolation_point(I.neighbour((dir_sons)i),level);
      }
    }
    
    // boundary points
    iterate_hash3 {
      I   = bo2point->ind;
      Iop = I.next(bo2point->direction,level);
      Add_interpolation_point(I,level);
      Add_interpolation_point(Iop,level);
    }
    
    // parallelization
    if(parallel_version) {
      Send_multi_grid_points_parallel(max_level-1);
    }

    Index3D my_lev_index;

    // second part: coarser points
    // -----------
    //    for(level=max_level-1;level>=1;--level) {

    for(level = max_level-1;
        level >= Give_my_coarsest_level() && level >1;
        --level) {
      my_lev_index = Give_my_level_index(level+1);
      
      iterate_hash1 {
        I = point1->Give_Index();
        if(I << my_lev_index) {
            typI = Give_type(I);

            if(typI>interior && I.Tiefe()==level) {
              Add_interpolation_point(I,level);
            }
            //      if(typI==multigrid && I.Tiefe()==level) {
            if(typI==multigrid && I.Tiefe()==level) {


              if(point1->finest_level>=level)
                Add_interpolation_point(I,level);
            }
          }
        }
      
      // parallelization
      if(parallel_version) {
        Send_multi_grid_points_parallel(level-1);
      }
    }
    
    // three part: make non-multigrid points exterior
    // -----------
    iterate_hash1 {
      I = point1->Give_Index();
      typI               = point1->Give_typ();
      if(typI==multigrid && I.Tiefe()<max_level && 
         point1->finest_level==0) {
        point1->typ=exterior;
      }
      point1->global_typ = point1->Give_typ();
    }
  }

  //Test_just_this();



  if(parallel_version) {
    // fourth part: find out responsibility
    //-------------
    int number_of_adjecent;
    Index3D adjecent_proc[8];
    int     adjecent_rank[8];

    int lowest_rank;
    bool found_responsibility;
    bool I_am_responsible;
    

    /*
     
     case i)

     *---*
     |  I|
     | + |
     |   |
     *---*

     case ii)

     *---*---*
     |   |   |
     | + I + |
     |   |   |
     *---*---*

     case iii)

     *---*---*
     |   |   |
     | + | + |
     |   |   |
     *---I---*
     |   |   |
     | + | + |
     |   |   |
     *---*---*

     +                 : is the center of the adjecent processors
     number_of_adjecent: is the number of thoses adjecent processors 
                         which exists


     */



    iterate_hash1 {
      // start
      I_am_responsible     = false;
      found_responsibility = false;
      
      I = point1->Give_Index();
      typI = Give_type(I);
      if(typI==multigrid) {
        // study this multigrid point I
        if(I<my_index) {
          // resonsibility case i)
          I_am_responsible     = true;
          found_responsibility = true;
        }
        else {
          // find processors which are adjecent to I
          number_of_adjecent = 0;
          /*
            for(p_iter=0;p_iter<hashtable_proc_leng;++p_iter) 
            for(point_proc=hashtable_proc_start[p_iter];point_proc!=NULL;
            point_proc=point_proc->next) {
          */
          iterate_hash_proc {
            Iproc=point_proc->Give_Index();
            if(Iproc.Tiefe()==n_parallel) {
              if(I<<Iproc) {
                // found one adjecent processor
                adjecent_proc[number_of_adjecent] = Iproc;
                adjecent_rank[number_of_adjecent] = 
                  point_proc->Give_num_proc();
                ++number_of_adjecent;
                if(I<Iproc) {
                  // resonsibility case i)
                  found_responsibility = true;
                }
              }
            }
            /*
              if(found_responsibility) {
              // stop iterate_hash_proc
              p_iter     = hashtable_proc_leng;
              point_proc = NULL;
              }
            */
          }

          if(developer_version) {
            if(found_responsibility==false &&
               I.Tiefe() >= n_parallel &&
               number_of_adjecent==0) {
              cout << " error 1 in Calc_multigrid_points_part2(); " << endl;
            }
          }

          if(found_responsibility == false) {
            if(I.Tiefe() >= n_parallel) {
              // responsibility case ii) 
              // point is contained in an edge or face of 
              // processor level n_parallel
              // (therefore not needed on coarser grids).
              
              // search lowest rank of adjecent processors
              lowest_rank = adjecent_rank[0];
              for(i=1;i<number_of_adjecent;++i) {
                if(lowest_rank > adjecent_rank[i])
                  lowest_rank = adjecent_rank[i];
              }
              found_responsibility = true;
              if(lowest_rank == my_rank) 
                I_am_responsible = true;
            }
            else {
              // resonsibility case iii)
              // point is corner of processor level n_parallel or coarser
              // (therefore also needed on coarser grids).
              if(developer_version) {
                if(Exists_proc_corner(I)==false) 
                  cout << " error 2 in Calc_multigrid_points_part2(); " 
                       << endl;
              }
              if(Give_proc_corner_number(I)==my_rank) {
                I_am_responsible = true;
              }
              else {
                I_am_responsible = false;
              }
            }
          }
        }
        if(I_am_responsible==false) {
          level = Give_finest_level(I);
          Set_point_typ(I,level,parallel_p);
        }
      }
    }
  }
}


void Grid_base::Add_interpolation_point(Index3D I, int level) {
  if(I.I_x().Tiefe()==level) {
    if(I.I_y().Tiefe()==level) {
      if(I.I_z().Tiefe()==level) {
        // middle point
        // edge between EST-WND
        Add_multigrid_point(I.next_EST(level),level-1);
        Add_multigrid_point(I.next_WND(level),level-1);
      }
      else {
        // TD-surface
        // edge between ES-WN
        Add_multigrid_point(I.next_ES(level),level-1);
        Add_multigrid_point(I.next_WN(level),level-1);
      }
    }
    else {
      if(I.I_z().Tiefe()==level) {
        // NS-surface
        // edge between WT-ED
        Add_multigrid_point(I.next_WT(level),level-1);
        Add_multigrid_point(I.next_ED(level),level-1);
      }
      else {
        // EW-edge
        Add_multigrid_point(I.next_E(level),level-1);
        Add_multigrid_point(I.next_W(level),level-1);
      }
    }
  }
  else {
    if(I.I_y().Tiefe()==level) {
      if(I.I_z().Tiefe()==level) {
        // EW-surface
        // edge between ST-ND
        Add_multigrid_point(I.next_ST(level),level-1);
        Add_multigrid_point(I.next_ND(level),level-1);
      }
      else {
        // NS-edge
        Add_multigrid_point(I.next_N(level),level-1);
        Add_multigrid_point(I.next_S(level),level-1);
      }
    }
    else {
      if(I.I_z().Tiefe()==level) {
        // TD-edge
        Add_multigrid_point(I.next_T(level),level-1);
        Add_multigrid_point(I.next_D(level),level-1);
      }
    }
  }
}


// schneller machen mit point1 anstatt I
void Grid_base::Add_multigrid_point(Index3D I, int level) {
  Pointtype typI;
  typI = Give_type(I);

  // hinzufuegen von Punkten die wegen der Restriktionsformel
  // zu grob sind
  if(typI==exterior) {
    Add_point(I,multigrid);
    typI = multigrid;
  }


  if(typI==multigrid) {
    Put_finest_level_minimal(I,level);
  }


  if(developer_version) {
    if(typI==exterior) {
      cout << " Error in Add_multigrid_point: my_rank: " << my_rank
           << " x: " << I.I_x().get()
           << " y: " << I.I_y().get()
           << " z: " << I.I_z().get()
           << " level: " << level
           << " Exists: " << Exists_Point(I) 
           << endl;
    }
  }
}




//  2. von alter version; muss mal erneuert werden





inline double Inverse_of_weight(int anz) {
  if(anz==0) return 0.0;
  else       return 1.0 / anz;
}


// cc: coarse cell
// fc: fine   cell
// co: corner
inline double Calc_trans_weight(dir_sons cc, dir_sons fc, dir_sons co,
                                double* weights) {
  int i,j,k;
  i = 2*  (cc&1)     +  (fc&1)     +  (co&1);
  j = 2* ((cc>>1)&1) + ((fc>>1)&1) + ((co>>1)&1);
  k = 2* ((cc>>2)&1) + ((fc>>2)&1) + ((co>>2)&1);

  if(i==0 || i==4 ||
     j==0 || j==4 ||
     k==0 || k==4) return 0.0;
  else 
    return weights[(i-1) + 3*(j-1) + 9*(k-1)];
}


void Grid::Calc_Near_bo_MG_Coefficients(Index3D I, int l, double weight) {
  list_bilding_nearb = new Point_mg_coeff_nearb(I,Give_variable(I,l+1),weight,
                                                list_bilding_nearb);

  for(int i=0;i<6;++i) {
    if(Exists_Bo_Point(I,(dir_3D)i)) {
      list_bilding_bo = new Point_mg_coeff_bo(I,(dir_3D)i,
                                              Give_variable(I,(dir_3D)i),
                                              weight,list_bilding_bo);
    }
  }
};

// Berechnung der Mehrgitterkoeffizienten
void Grid::Calc_MG_Coefficients() {
  P_nearb    *it_n;
  Index3D I, next, next_cell;
  int color;



  cout << " Error! Do not use Calc_MG_Coefficients()!" << endl;

  for(int lev=Max_level()-1;lev>=Min_level();--lev) {
   for(color=0;color<8;++color)
    for(it_n=Start_P_nearb(lev,color);it_n!=NULL;it_n=it_n->Next()) {
    I = it_n->Ind();

    // Start:
    list_bilding_nearb = NULL;
    list_bilding_bo = NULL;

    // M:
    list_bilding_nearb = new Point_mg_coeff_nearb(I,Give_variable(I,lev+1),
                                                  1.0,
                                                  list_bilding_nearb);
    weight_27[trans27(MOrt,MOrt,MOrt)] = 1.0;

    /*
      BoCeData* bocedata;
      int k;
      int i,j;
      int Inv_weight;
      double weight;
      double* UU;
    // edges:
    // N 
    next = I.next_N(lev+1);
    Inv_weight = Give_Prolongation_weight(next);
    weight = Inverse_of_weight(Inv_weight);
    if(Inv_weight!=0) 
      Calc_Near_bo_MG_Coefficients(next,lev,weight);
    weight_27[trans27(MOrt,ROrt,MOrt)] = weight;
    // S
    next = I.next_S(lev+1);
    Inv_weight = Give_Prolongation_weight(next);
    weight = Inverse_of_weight(Inv_weight);
    if(Inv_weight!=0)   
      Calc_Near_bo_MG_Coefficients(next,lev,weight);
    weight_27[trans27(MOrt,LOrt,MOrt)] = weight;
    // E
    next = I.next_E(lev+1);
    Inv_weight = Give_Prolongation_weight(next);
    weight = Inverse_of_weight(Inv_weight);
    if(Inv_weight!=0)   
      Calc_Near_bo_MG_Coefficients(next,lev,weight);
    weight_27[trans27(ROrt,MOrt,MOrt)] = weight;
    // W
    next = I.next_W(lev+1);
    Inv_weight = Give_Prolongation_weight(next);
    weight = Inverse_of_weight(Inv_weight);
    if(Inv_weight!=0)   
      Calc_Near_bo_MG_Coefficients(next,lev,weight);
    weight_27[trans27(LOrt,MOrt,MOrt)] = weight;
    // T
    next = I.next_T(lev+1);
    Inv_weight = Give_Prolongation_weight(next);
    weight = Inverse_of_weight(Inv_weight);
    if(Inv_weight!=0)   
      Calc_Near_bo_MG_Coefficients(next,lev,weight);
    weight_27[trans27(MOrt,MOrt,ROrt)] = weight;
    // D
    next = I.next_D(lev+1);
    Inv_weight = Give_Prolongation_weight(next);
    weight = Inverse_of_weight(Inv_weight);
    if(Inv_weight!=0)   
      Calc_Near_bo_MG_Coefficients(next,lev,weight);
    weight_27[trans27(MOrt,MOrt,LOrt)] = weight;


    // classic faces:
    // ND
    next = I.next_ND(lev+1);
    Inv_weight = Give_Prolongation_weight(next);
    weight = Inverse_of_weight(Inv_weight);
    if(Inv_weight==4)
      list_bilding_nearb = 
        new Point_mg_coeff_nearb(next,Give_variable(next,lev+1),0.5,
                                 list_bilding_nearb);
    if(Inv_weight!=0)   
      Calc_Near_bo_MG_Coefficients(next,lev,weight);
    weight_27[trans27(MOrt,ROrt,LOrt)] = weight;
    // WN
    next = I.next_WN(lev+1);
    Inv_weight = Give_Prolongation_weight(next);
    weight = Inverse_of_weight(Inv_weight);
    if(Inv_weight==4)
      list_bilding_nearb = 
        new Point_mg_coeff_nearb(next,Give_variable(next,lev+1),0.5,
                                 list_bilding_nearb);
    if(Inv_weight!=0)   
      Calc_Near_bo_MG_Coefficients(next,lev,weight);
    weight_27[trans27(LOrt,ROrt,MOrt)] = weight;
    // WT
    next = I.next_WT(lev+1);
    Inv_weight = Give_Prolongation_weight(next);
    weight = Inverse_of_weight(Inv_weight);
    if(Inv_weight==4)
      list_bilding_nearb = 
        new Point_mg_coeff_nearb(next,Give_variable(next,lev+1),0.5,
                                 list_bilding_nearb);
    if(Inv_weight!=0)   
      Calc_Near_bo_MG_Coefficients(next,lev,weight);
    weight_27[trans27(LOrt,MOrt,ROrt)] = weight;
    // ED
    next = I.next_ED(lev+1);
    Inv_weight = Give_Prolongation_weight(next);
    weight = Inverse_of_weight(Inv_weight);
    if(Inv_weight==4)
      list_bilding_nearb = 
        new Point_mg_coeff_nearb(next,Give_variable(next,lev+1),0.5,
                                 list_bilding_nearb);
    if(Inv_weight!=0)   
      Calc_Near_bo_MG_Coefficients(next,lev,weight);
    weight_27[trans27(ROrt,MOrt,LOrt)] = weight;
    // ST
    next = I.next_ST(lev+1);
    Inv_weight = Give_Prolongation_weight(next);
    weight = Inverse_of_weight(Inv_weight);
    if(Inv_weight==4)
      list_bilding_nearb = 
        new Point_mg_coeff_nearb(next,Give_variable(next,lev+1),0.5,
                                 list_bilding_nearb);
    if(Inv_weight!=0)   
      Calc_Near_bo_MG_Coefficients(next,lev,weight);
    weight_27[trans27(MOrt,LOrt,ROrt)] = weight;
    // ES
    next = I.next_ES(lev+1);
    Inv_weight = Give_Prolongation_weight(next);
    weight = Inverse_of_weight(Inv_weight);
    if(Inv_weight==4)
      list_bilding_nearb = 
        new Point_mg_coeff_nearb(next,Give_variable(next,lev+1),0.5,
                                 list_bilding_nearb);
    if(Inv_weight!=0)   
      Calc_Near_bo_MG_Coefficients(next,lev,weight);
    weight_27[trans27(ROrt,LOrt,MOrt)] = weight;



    // only near boundary faces:
    // SD
    next = I.next_SD(lev+1);
    Inv_weight = Give_Prolongation_weight(next);
    weight = Inverse_of_weight(Inv_weight);
    if(Inv_weight!=0 && Inv_weight!=4)   
      Calc_Near_bo_MG_Coefficients(next,lev,weight);
    weight_27[trans27(MOrt,LOrt,LOrt)] = weight;
    // NT
    next = I.next_NT(lev+1);
    Inv_weight = Give_Prolongation_weight(next);
    weight = Inverse_of_weight(Inv_weight);
    if(Inv_weight!=0 && Inv_weight!=4)   
      Calc_Near_bo_MG_Coefficients(next,lev,weight);
    weight_27[trans27(MOrt,ROrt,ROrt)] = weight;
    // EN
    next = I.next_EN(lev+1);
    Inv_weight = Give_Prolongation_weight(next);
    weight = Inverse_of_weight(Inv_weight);
    if(Inv_weight!=0 && Inv_weight!=4)   
      Calc_Near_bo_MG_Coefficients(next,lev,weight);
    weight_27[trans27(ROrt,ROrt,MOrt)] = weight;
    // WS
    next = I.next_WS(lev+1);
    Inv_weight = Give_Prolongation_weight(next);
    weight = Inverse_of_weight(Inv_weight);
    if(Inv_weight!=0 && Inv_weight!=4)   
      Calc_Near_bo_MG_Coefficients(next,lev,weight);
    weight_27[trans27(LOrt,LOrt,MOrt)] = weight;
    // ET
    next = I.next_ET(lev+1);
    Inv_weight = Give_Prolongation_weight(next);
    weight = Inverse_of_weight(Inv_weight);
    if(Inv_weight!=0 && Inv_weight!=4)   
      Calc_Near_bo_MG_Coefficients(next,lev,weight);
    weight_27[trans27(ROrt,MOrt,ROrt)] = weight;
    // WD
    next = I.next_WD(lev+1);
    Inv_weight = Give_Prolongation_weight(next);
    weight = Inverse_of_weight(Inv_weight);
    if(Inv_weight!=0 && Inv_weight!=4)   
      Calc_Near_bo_MG_Coefficients(next,lev,weight);
    weight_27[trans27(LOrt,MOrt,LOrt)] = weight;


    // classic cells:
    // EST
    next = I.next_EST(lev+1);
    Inv_weight = Give_Prolongation_weight(next);
    weight = Inverse_of_weight(Inv_weight);
    if(Inv_weight==8)
      list_bilding_nearb = 
        new Point_mg_coeff_nearb(next,Give_variable(next,lev+1),0.5,
                                 list_bilding_nearb);
    if(Inv_weight!=0)   
      Calc_Near_bo_MG_Coefficients(next,lev,weight);
    weight_27[trans27(ROrt,LOrt,ROrt)] = weight;

    // WND
    next = I.next_WND(lev+1);
    Inv_weight = Give_Prolongation_weight(next);
    weight = Inverse_of_weight(Inv_weight);
    if(Inv_weight==8)
      list_bilding_nearb = 
        new Point_mg_coeff_nearb(next,Give_variable(next,lev+1),0.5,
                                 list_bilding_nearb);
    if(Inv_weight!=0)   
      Calc_Near_bo_MG_Coefficients(next,lev,weight);
    weight_27[trans27(LOrt,ROrt,LOrt)] = weight;

    // only near boundary cells:
    // ESD
    next = I.next_ESD(lev+1);
    Inv_weight = Give_Prolongation_weight(next);
    weight = Inverse_of_weight(Inv_weight);
    if(Inv_weight!=0 && Inv_weight!=8)   
      Calc_Near_bo_MG_Coefficients(next,lev,weight);
    weight_27[trans27(ROrt,LOrt,LOrt)] = weight;
    // ENT
    next = I.next_ENT(lev+1);
    Inv_weight = Give_Prolongation_weight(next);
    weight = Inverse_of_weight(Inv_weight);
    if(Inv_weight!=0 && Inv_weight!=8)   
      Calc_Near_bo_MG_Coefficients(next,lev,weight);
    weight_27[trans27(ROrt,ROrt,ROrt)] = weight;
    // END
    next = I.next_END(lev+1);
    Inv_weight = Give_Prolongation_weight(next);
    weight = Inverse_of_weight(Inv_weight);
    if(Inv_weight!=0 && Inv_weight!=8)   
      Calc_Near_bo_MG_Coefficients(next,lev,weight);
    weight_27[trans27(ROrt,ROrt,LOrt)] = weight;
    // WSD
    next = I.next_WSD(lev+1);
    Inv_weight = Give_Prolongation_weight(next);
    weight = Inverse_of_weight(Inv_weight);
    if(Inv_weight!=0 && Inv_weight!=8)   
      Calc_Near_bo_MG_Coefficients(next,lev,weight);
    weight_27[trans27(LOrt,LOrt,LOrt)] = weight;
    // WNT
    next = I.next_WNT(lev+1);
    Inv_weight = Give_Prolongation_weight(next);
    weight = Inverse_of_weight(Inv_weight);
    if(Inv_weight!=0 && Inv_weight!=8)   
      Calc_Near_bo_MG_Coefficients(next,lev,weight);
    weight_27[trans27(LOrt,ROrt,ROrt)] = weight;
    // WST
    next = I.next_WST(lev+1);
    Inv_weight = Give_Prolongation_weight(next);
    weight = Inverse_of_weight(Inv_weight);
    if(Inv_weight!=0 && Inv_weight!=8)   
      Calc_Near_bo_MG_Coefficients(next,lev,weight);
    weight_27[trans27(LOrt,LOrt,ROrt)] = weight;

    // bocell point (Randzellfreiheitsgrad)
    for(i=0;i<8;++i) {    // go to middle point of neighbour coarse cells
      next_cell = I.next((dir_sons)i,lev+1); 
      for(j=0;j<8;++j) {    // go to middle point of neighbour fine cells
        next = next_cell.next((dir_sons)i,lev+2); 
        bocedata = Give_Bo_cell(next);
        if(bocedata!=NULL) {
          if(bocedata->Exists_bocellpoint()) {
            //    cout << " fine_bo_cell " << endl;
            UU = Give_vector_ablage();
            
            for(k=0;k<8;++k) {
              UU[k] = Calc_trans_weight((dir_sons)i,(dir_sons)j,(dir_sons)k,
                                        weight_27);
            }
            list_bilding_nearb = new Point_mg_coeff_nearb(I,
                      bocedata->vars[bocedata->Give_number_points()],
                      bocedata->Weight_bocellpoint() *
                        bocedata->Interpolation_interior_corners(UU),
                      list_bilding_nearb);
          }
        }
      }
    }
    it_n->Set_first_mg_coeff_nearb(list_bilding_nearb);
    it_n->Set_first_mg_coeff_bo(list_bilding_bo);
    */

    }
  }
};

---