src/expde/extemp/locsum.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
//

// ------------------------------------------------------------
// locsum.cc
//
// ------------------------------------------------------------

// Include level 0:
#ifdef COMP_GNUOLD
 #include <iostream.h>
 #include <fstream.h>
 #include <time.h>
#else
 #include <iostream>
 #include <fstream>
 #include <ctime>
#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 "variable.h"
#include "locsum.h"


#define sum_up_27_test(ppoi) ppoi->varM(grid,level)[num_sum_var];

#define sum_up_27(ppoi) ppoi->varM(grid,level)[num_sum_var] + \
                        ppoi->varEND(grid,level)[num_sum_var] + \
                        ppoi->varWND(grid,level)[num_sum_var] + \
                        ppoi->varESD(grid,level)[num_sum_var] + \
                        ppoi->varWSD(grid,level)[num_sum_var] + \
                        ppoi->varENT(grid,level)[num_sum_var] + \
                        ppoi->varWNT(grid,level)[num_sum_var] + \
                        ppoi->varEST(grid,level)[num_sum_var] + \
                        ppoi->varWST(grid,level)[num_sum_var] + \
                        ppoi->varND(grid,level)[num_sum_var] + \
                        ppoi->varSD(grid,level)[num_sum_var] + \
                        ppoi->varNT(grid,level)[num_sum_var] + \
                        ppoi->varST(grid,level)[num_sum_var] + \
                        ppoi->varED(grid,level)[num_sum_var] + \
                        ppoi->varWD(grid,level)[num_sum_var] + \
                        ppoi->varET(grid,level)[num_sum_var] + \
                        ppoi->varWT(grid,level)[num_sum_var] + \
                        ppoi->varEN(grid,level)[num_sum_var] + \
                        ppoi->varWN(grid,level)[num_sum_var] + \
                        ppoi->varES(grid,level)[num_sum_var] + \
                        ppoi->varWS(grid,level)[num_sum_var] + \
                        ppoi->varN(grid,level)[num_sum_var] + \
                        ppoi->varS(grid,level)[num_sum_var] + \
                        ppoi->varE(grid,level)[num_sum_var] + \
                        ppoi->varW(grid,level)[num_sum_var] + \
                        ppoi->varT(grid,level)[num_sum_var] + \
                        ppoi->varD(grid,level)[num_sum_var];

#define set_next_27(ppoi) ppoi->varEND(grid,level)[number_variable] = sum; \
                          ppoi->varWND(grid,level)[number_variable] = sum; \
                          ppoi->varESD(grid,level)[number_variable] = sum; \
                          ppoi->varWSD(grid,level)[number_variable] = sum; \
                          ppoi->varENT(grid,level)[number_variable] = sum; \
                          ppoi->varWNT(grid,level)[number_variable] = sum; \
                          ppoi->varEST(grid,level)[number_variable] = sum; \
                          ppoi->varWST(grid,level)[number_variable] = sum; \
                          ppoi->varND(grid,level)[number_variable] = sum; \
                          ppoi->varSD(grid,level)[number_variable] = sum; \
                          ppoi->varNT(grid,level)[number_variable] = sum; \
                          ppoi->varST(grid,level)[number_variable] = sum; \
                          ppoi->varED(grid,level)[number_variable] = sum; \
                          ppoi->varWD(grid,level)[number_variable] = sum; \
                          ppoi->varET(grid,level)[number_variable] = sum; \
                          ppoi->varWT(grid,level)[number_variable] = sum; \
                          ppoi->varEN(grid,level)[number_variable] = sum; \
                          ppoi->varWN(grid,level)[number_variable] = sum; \
                          ppoi->varES(grid,level)[number_variable] = sum; \
                          ppoi->varWS(grid,level)[number_variable] = sum; \
                          ppoi->varN(grid,level)[number_variable] = sum; \
                          ppoi->varS(grid,level)[number_variable] = sum; \
                          ppoi->varE(grid,level)[number_variable] = sum; \
                          ppoi->varW(grid,level)[number_variable] = sum; \
                          ppoi->varT(grid,level)[number_variable] = sum; \
                          ppoi->varD(grid,level)[number_variable] = sum;


// some first functions
LocalSumObj Local_Sum(int i, int j, int k, Variable& v) {
  return LocalSumObj(i,j,k, v.Number_variable(),v.Level()); 
}




void Variable::operator<<= (const LocalSumObj& local_sum_object) {
  int I,J,K;
  int p,q,r;
  dir_sons color;
  int level;
  int num_sum_var;
  double sum;
  Index3D ind;

  // 1. Step: get values
  I = local_sum_object.i();
  J = local_sum_object.j();
  K = local_sum_object.k();

  num_sum_var = local_sum_object.Number_variable();
  level = local_sum_object.Level();

  // 2. Step: test error in Parameter
  if(I<0 || I>3 ||
     J<0 || J>3 ||
     K<0 || K>3) {
    cout << " Parameter Error in Local_Sum! " << endl;
  }
  else {
    // 3. Step calculate color
    color = (dir_sons)((I&1) + 2 * (J&1) + 4 * (K&1));
    color = reordering_formula(color);

    /*
    cout << " jjjjj: " << color 
               << " I: " << (I&1)
               << " J: " << (J&1)
               << " K: " << (K&1) 
               << " level: " << level << endl;
    */

    if(level=Max_Level()) { // !!! on finest grid
      // 4. Step: Sum up
      // innere Punkte
      for(iter_i=grid->Start_P_interior(level,color);
          iter_i!=NULL;iter_i=iter_i->Next()) {
        ind = iter_i->Ind();
        if(ind.Is_color_four(I,J,K,level)) {
          
          sum = sum_up_27(iter_i);
          for(p=0;p<3;++p) for(q=0;q<3;++q) for(r=0;r<3;++r) {
            sum = sum + grid->sum_rest_at_point(ind.next((Ort1D)p,
                                                         (Ort1D)q,
                                                         (Ort1D)r,level),
                                                num_sum_var);
          }
          sum = sum + grid->sum_of_cellp(ind,num_sum_var);

          iter_i->varM(grid,level)[number_variable] = sum;      
          
        }
      }
      
      // randnahe Punkte
      for(iter_n=grid->Start_P_nearb(level,color);
          iter_n!=NULL;iter_n=iter_n->Next()) {
        ind = iter_n->Ind();
        if(ind.Is_color_four(I,J,K,level)) {
          sum = 0.0;
          for(p=0;p<3;++p) for(q=0;q<3;++q) for(r=0;r<3;++r) {
            sum = sum + grid->sum_at_point(ind.next((Ort1D)p,
                                                    (Ort1D)q,
                                                    (Ort1D)r,level),
                                           num_sum_var,level);
          }
          sum = sum + grid->sum_of_cellp(ind,num_sum_var);
          //    sum = sum_up_27_test(iter_n);
          
          iter_n->varM(grid,level)[number_variable] = sum;
        }
      }
      // 5. Step: Set next points
      // innere Punkte
      for(iter_i=grid->Start_P_interior(level,color);
          iter_i!=NULL;iter_i=iter_i->Next()) {
        ind = iter_i->Ind();
        if(ind.Is_color_four(I,J,K,level)) {
          sum = iter_i->varM(grid,level)[number_variable];      
          set_next_27(iter_i);
          for(p=0;p<3;++p) for(q=0;q<3;++q) for(r=0;r<3;++r) {
            grid->set_rest_at_point(ind.next((Ort1D)p,
                                             (Ort1D)q,
                                             (Ort1D)r,level),
                                    number_variable,sum);
            grid->set_at_cellp(ind,number_variable,sum);
          }
        }
      }
      
      // randnahe Punkte
      for(iter_n=grid->Start_P_nearb(level,color);
          iter_n!=NULL;iter_n=iter_n->Next()) {
        ind = iter_n->Ind();
        if(ind.Is_color_four(I,J,K,level)) {
          sum = iter_n->varM(grid,level)[number_variable];      
          for(p=0;p<3;++p) for(q=0;q<3;++q) for(r=0;r<3;++r) {
            grid->set_at_point(ind.next((Ort1D)p,
                                        (Ort1D)q,
                                        (Ort1D)r,level),
                               number_variable,sum,level);
            grid->set_at_cellp(ind,number_variable,sum);
          }
        }
      }
    }
    else { // !!! on coarser grids
      // 4. Step: Sum up
      // innere Punkte
      for(iter_i=grid->Start_P_interior(level,color);
          iter_i!=NULL;iter_i=iter_i->Next()) {
        ind = iter_i->Ind();
        if(ind.Is_color_four(I,J,K,level)) {      
          sum = sum_up_27(iter_i);
          iter_i->varM(grid,level)[number_variable] = sum;                
        }
      }
      
      // randnahe Punkte
      for(iter_n=grid->Start_P_nearb(level,color);
          iter_n!=NULL;iter_n=iter_n->Next()) {
        ind = iter_n->Ind();
        if(ind.Is_color_four(I,J,K,level)) {
          sum = 0.0;
          for(p=0;p<3;++p) for(q=0;q<3;++q) for(r=0;r<3;++r) {
            sum = sum + grid->sum_at_point(ind.next((Ort1D)p,
                                                    (Ort1D)q,
                                                    (Ort1D)r,level),
                                           num_sum_var,level);
          }
          iter_n->varM(grid,level)[number_variable] = sum;
        }
      }

      // exterior grobe Randpunkte
      for(iter_n=grid->Start_P_exteri(level,color);
          iter_n!=NULL;iter_n=iter_n->Next()) {
        ind = iter_n->Ind();
        if(ind.Is_color_four(I,J,K,level)) {
          sum = 0.0;
          for(p=0;p<3;++p) for(q=0;q<3;++q) for(r=0;r<3;++r) {
            sum = sum + grid->sum_at_point(ind.next((Ort1D)p,
                                                    (Ort1D)q,
                                                    (Ort1D)r,level),
                                           num_sum_var,level);
          }
          iter_n->varM(grid,level)[number_variable] = sum;
        }
      }
      // 5. Step: Sum up
      // innere Punkte
      for(iter_i=grid->Start_P_interior(level,color);
          iter_i!=NULL;iter_i=iter_i->Next()) {
        ind = iter_i->Ind();
        if(ind.Is_color_four(I,J,K,level)) {
          sum = iter_i->varM(grid,level)[number_variable];
          set_next_27(iter_i);     
        }
      }
      
      // randnahe Punkte
      for(iter_n=grid->Start_P_nearb(level,color);
          iter_n!=NULL;iter_n=iter_n->Next()) {
        ind = iter_n->Ind();
        if(ind.Is_color_four(I,J,K,level)) {
          sum = iter_n->varM(grid,level)[number_variable];
          for(p=0;p<3;++p) for(q=0;q<3;++q) for(r=0;r<3;++r) {
            grid->set_at_point(ind.next((Ort1D)p,
                                        (Ort1D)q,
                                        (Ort1D)r,level),
                               number_variable,sum,level);
          }                                                                
        }
      }

      // exterior grobe Randpunkte
      for(iter_n=grid->Start_P_exteri(level,color);
          iter_n!=NULL;iter_n=iter_n->Next()) {
        ind = iter_n->Ind();
        if(ind.Is_color_four(I,J,K,level)) {
          sum = iter_n->varM(grid,level)[number_variable];
          for(p=0;p<3;++p) for(q=0;q<3;++q) for(r=0;r<3;++r) {
            grid->set_at_point(ind.next((Ort1D)p,
                                        (Ort1D)q,
                                        (Ort1D)r,level),
                               number_variable,sum,level);
          }
        }
      }
    }

    
  }
}


// member functions 

LocalSumObj::LocalSumObj(int i, int j, int k, int num_var, int lev) {
  I=i;
  J=j;
  K=k;

  number_var = num_var;
  level = lev;
}


double Grid_base::sum_rest_at_point(Index3D ind, int num_sum_variable) {
  double sum;
  sum = 0.0;
  if(Give_type(ind)>interior) {
    for(int d=0;d<6;d++) {
      if(Exists_Bo_Point(ind,(dir_3D)d)) {
        sum = sum + Give_variable(ind,(dir_3D)d)[num_sum_variable];
      }
    }
  }
  return sum;
}

double Grid_base::sum_at_point(Index3D ind, int num_sum_variable, int level) {
  double sum;

  if(Give_type(ind)>=interior) {
    sum = Give_variable(ind,level)[num_sum_variable];
    if(level == max_level) {
      for(int d=0;d<6;d++) {
        if(Exists_Bo_Point(ind,(dir_3D)d)) {
          sum = sum + Give_variable(ind,(dir_3D)d)[num_sum_variable];
        }
      }
    }
  }
  else {
    sum = 0.0;
  }
  return sum;
}

double Grid_base::sum_of_cellp(Index3D ind, int num_sum_variable) {
  double sum;
  Index3D In;
  Index3D I;
  double* vars;

  sum = 0.0;
  for(int i=0;i<8;++i) { // go to next corners
    In = ind.next((dir_sons)i,max_level);
    for(int j=0;j<8;++j) { // go to next cells
     I = In.next((dir_sons)j,max_level+1);
     vars = Give_variable_cellpoi_slow(I);
     if(vars!=NULL) {
       sum = sum + vars[num_sum_variable];
     }
    }
  }
  return sum;
}


void Grid_base::set_rest_at_point(Index3D ind, int num_variable, 
                                  double sum) {
  if(Give_type(ind)>interior) {
    for(int d=0;d<6;d++) {
      if(Exists_Bo_Point(ind,(dir_3D)d)) {
        Give_variable(ind,(dir_3D)d)[num_variable] = sum;
      }
    }
  }
}

void Grid_base::set_at_point(Index3D ind, int num_variable,
                             double sum, int level) {
  if(Give_type(ind)>=interior) {
    Give_variable(ind,level)[num_variable] = sum;
    if(level == max_level) {
      for(int d=0;d<6;d++) {
        if(Exists_Bo_Point(ind,(dir_3D)d)) {
          Give_variable(ind,(dir_3D)d)[num_variable] = sum;
        }
      }
    }
  }
}

void Grid_base::set_at_cellp(Index3D ind, int number_variable, double sum) {
  Index3D In;
  Index3D I;
  double* vars;

  for(int i=0;i<8;++i) { // go to next corners
    In = ind.next((dir_sons)i,max_level);
    for(int j=0;j<8;++j) { // go to next cells
     I = In.next((dir_sons)j,max_level+1);
     vars = Give_variable_cellpoi_slow(I);
     if(vars!=NULL) {
       vars[number_variable] = sum;
     }
    }
  }
}

---