src/expde/formulas/res_sten.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
//

// ------------------------------------------------------------
// res_sten.cc
//
// ------------------------------------------------------------

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


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

// Include level special
#include "../formulas/loc_sten.h"


// This file contains the restriction formula for
// local 64- Stiffness matrices
// 1. formula for interior points
// 2. formula for boundary points


// 1. formula for interior points


inline int number_coarse(dir_sons i) {
  return trans_cellP_27(i,i);
}

void Set_coarse_function(int* U, dir_sons v) {
  int k, trans, d;
  dir_sons next;


  // part 0: set zero
  //--------
  for(k=0;k<27;++k) U[k]=0;

  // part 1: corners of coarse cell
  //--------
  for(k=0;k<8;++k) { // iteration over cells
    if(k==v) U[number_coarse((dir_sons)k)]=2;
    else     U[number_coarse((dir_sons)k)]=0;
  }

  // part 2: middle of coarse cell
  //--------
  U[trans27(MOrt,MOrt,MOrt)] = (U[number_coarse(ESTd)] +
                                U[number_coarse(WNDd)]) / 2;

  // part 3: faces of coarse cell
  //--------
  // TD-surface
  // edge between ES-WN
  U[trans27(MOrt,MOrt,LOrt)] = (U[number_coarse(ESDd)] +
                                U[number_coarse(WNDd)]) / 2;
  U[trans27(MOrt,MOrt,ROrt)] = (U[number_coarse(ESTd)] +
                                U[number_coarse(WNTd)]) / 2;
  // NS-surface
  // edge between WT-ED
  U[trans27(MOrt,LOrt,MOrt)] = (U[number_coarse(WSTd)] +
                                U[number_coarse(ESDd)]) / 2;
  U[trans27(MOrt,ROrt,MOrt)] = (U[number_coarse(WNTd)] +
                                U[number_coarse(ENDd)]) / 2;
  // EW-surface
  // edge between ST-ND
  U[trans27(LOrt,MOrt,MOrt)] = (U[number_coarse(WSTd)] +
                                U[number_coarse(WNDd)]) / 2;
  U[trans27(ROrt,MOrt,MOrt)] = (U[number_coarse(ESTd)] +
                                U[number_coarse(ENDd)]) / 2;

  // part 4: edges of coarse cell
  //--------
  for(k=0;k<8;++k) { // iteration over cells
    for(d=0;d<3;++d) {  // iteration over directions
      next = Next_corner((dir_sons)k,(dir_3D)(2*d));
      trans = trans_cellP_27((dir_sons)k,next);
      U[trans] = U[trans] + U[number_coarse((dir_sons)k)] / 2;
    }
  }
}


void Restrict_local_stiffness_matrix(double*  Stencil_coarse,
                                     double** Stencils_fine,
                                     int* coarse_grid_functions) {

  int i,j,k,l,m;
  int* v_function;
  int* u_function;
  int u_f, v_f;

  // 1.Part: Set zero
  for(i=0;i<8;++i) for(j=0;j<8;++j) Stencil_coarse[j+8*i] = 0.0;
  
  // 2.Part: Iterate
  for(i=0;i<8;++i) { // iteration for U
    u_function = &(coarse_grid_functions[i*27]);
    for(j=0;j<8;++j) { // iteration for V
      v_function = &(coarse_grid_functions[j*27]);
      for(k=0;k<8;++k) { // iteration over cells
        for(l=0;l<8;++l) { // iteration over corners of cell for U
          u_f = u_function[trans_cellP_27((dir_sons)k,(dir_sons)l)];
          if(u_f != 0) for(m=0;m<8;++m) { // iteration over corners of cell for V
            v_f = v_function[trans_cellP_27((dir_sons)k,(dir_sons)m)];
            if(v_f!=0)
              Stencil_coarse[j+8*i] = Stencil_coarse[j+8*i] +
                (double)u_f * (double)v_f * 0.25 *
                Stencils_fine[k][m+8*l];
          }
        }
      }
    }
  }

}




// 2. formula for boundary points


// cid: corner in domain
void Set_coarse_function(double U[27], dir_sons v, bool* cid) {
  int i,j,k;
  dir_sons corner;
  int anz;
  double sum;

  // part 0: set zero
  //--------
  for(k=0;k<27;++k) U[k]=0.0;

  // part 1: corners of coarse cell
  //--------
  for(k=0;k<8;++k) { // iteration over cells
    if(k==v) U[number_coarse((dir_sons)k)]=1.0;
    else     U[number_coarse((dir_sons)k)]=0.0;
  }

  // part 2: middle of coarse cell
  //--------
  anz = 0;
  for(i=0;i<2;++i) for(j=0;j<2;++j) for(k=0;k<2;++k) {
    corner = dir_sons_from_three_LR((dir1D)i,(dir1D)j,(dir1D)k);
    if(cid[corner]) {
      ++anz;
      U[trans27(MOrt,MOrt,MOrt)] +=
        U[number_coarse(corner)];
    }
  }
  if(anz!=0) U[trans27(MOrt,MOrt,MOrt)] = U[trans27(MOrt,MOrt,MOrt)] / anz;
  else       U[trans27(MOrt,MOrt,MOrt)] = 0.0;

  // part 3: faces of coarse cell
  //--------
  // part 3.1: interpolation for a interior face
  //----------
  // TD-surface
  // edge between ES-WN
  U[trans27(MOrt,MOrt,LOrt)] = 0.5*(U[number_coarse(ESDd)] +
                                    U[number_coarse(WNDd)]);
  U[trans27(MOrt,MOrt,ROrt)] = 0.5*(U[number_coarse(ESTd)] +
                                    U[number_coarse(WNTd)]);
  // NS-surface
  // edge between WT-ED
  U[trans27(MOrt,LOrt,MOrt)] = 0.5*(U[number_coarse(WSTd)] +
                                    U[number_coarse(ESDd)]);
  U[trans27(MOrt,ROrt,MOrt)] = 0.5*(U[number_coarse(WNTd)] +
                                    U[number_coarse(ENDd)]);
  // EW-surface
  // edge between ST-ND
  U[trans27(LOrt,MOrt,MOrt)] = 0.5*(U[number_coarse(WSTd)] +
                                    U[number_coarse(WNDd)]);
  U[trans27(ROrt,MOrt,MOrt)] = 0.5*(U[number_coarse(ESTd)] +
                                    U[number_coarse(ENDd)]);
  // part 3.2: interpolation at a face near the boundary
  // TD-surface
  for(k=0;k<2;++k) {
    anz=0;
    sum=0.0;
    for(i=0;i<2;++i) for(j=0;j<2;++j) {
      corner = dir_sons_from_three_LR((dir1D)i,(dir1D)j,(dir1D)k);
      if(cid[corner]) {
        ++anz;
        sum += U[number_coarse(corner)];
      }
    }
    if(anz!=0) {
      if(anz!=4) // overwrite 3.1
        U[trans27(MOrt,MOrt,(Ort1D)(2*k))] = sum / anz;
    }
    else U[trans27(MOrt,MOrt,(Ort1D)(2*k))] = 0.0;
  }
  // NS-surface
  for(j=0;j<2;++j) {
    anz=0;
    sum=0.0;
    for(i=0;i<2;++i) for(k=0;k<2;++k) {
      corner = dir_sons_from_three_LR((dir1D)i,(dir1D)j,(dir1D)k);
      if(cid[corner]) {
        ++anz;
        sum += U[number_coarse(corner)];
      }
    }
    if(anz!=0) {
      if(anz!=4)  // overwrite 3.1
        U[trans27(MOrt,(Ort1D)(2*j),MOrt)] = sum / anz;
    }
    else U[trans27(MOrt,(Ort1D)(2*j),MOrt)] = 0.0;
  }
  // EW-surface
  for(i=0;i<2;++i) {
    anz=0;
    sum=0.0;
    for(j=0;j<2;++j) for(k=0;k<2;++k) {
      corner = dir_sons_from_three_LR((dir1D)i,(dir1D)j,(dir1D)k);
      if(cid[corner]) {
        ++anz;
        sum += U[number_coarse(corner)];
      }
    }
    if(anz!=0) {
      if(anz!=4)  // overwrite 3.1
        U[trans27((Ort1D)(2*i),MOrt,MOrt)] = sum / anz;
    }
    else U[trans27((Ort1D)(2*i),MOrt,MOrt)] = 0.0;
  }

  // part 4: edges of coarse cell
  //--------
  // EW-edges
  for(j=0;j<2;++j) for(k=0;k<2;++k) {
    anz=0;
    for(i=0;i<2;++i) {
      corner = dir_sons_from_three_LR((dir1D)i,(dir1D)j,(dir1D)k);
      if(cid[corner]) {
        ++anz;
        U[trans27(MOrt,(Ort1D)(2*j),(Ort1D)(2*k))] +=
          U[number_coarse(corner)];
      }
    }
    if(anz!=0) 
      U[trans27(MOrt,(Ort1D)(2*j),(Ort1D)(2*k))] = 
        U[trans27(MOrt,(Ort1D)(2*j),(Ort1D)(2*k))] / anz;
    else U[trans27(MOrt,(Ort1D)(2*j),(Ort1D)(2*k))] = 0.0;
  }

  // NS-edges
  for(i=0;i<2;++i) for(k=0;k<2;++k) {
    anz=0;
    for(j=0;j<2;++j) {
      corner = dir_sons_from_three_LR((dir1D)i,(dir1D)j,(dir1D)k);
      if(cid[corner]) {
        ++anz;
        U[trans27((Ort1D)(2*i),MOrt,(Ort1D)(2*k))] +=
          U[number_coarse(corner)];
      }
    }
    if(anz!=0) {
      U[trans27((Ort1D)(2*i),MOrt,(Ort1D)(2*k))] = 
        U[trans27((Ort1D)(2*i),MOrt,(Ort1D)(2*k))] / anz;
    }
    else U[trans27((Ort1D)(2*i),MOrt,(Ort1D)(2*k))] = 0.0;
  }

  // TD-edges
  for(i=0;i<2;++i) for(j=0;j<2;++j) {
    anz=0;
    for(k=0;k<2;++k) {
      corner = dir_sons_from_three_LR((dir1D)i,(dir1D)j,(dir1D)k);
      if(cid[corner]) {
        ++anz;
        U[trans27((Ort1D)(2*i),(Ort1D)(2*j),MOrt)] +=
          U[number_coarse(corner)];
      }
    }
    if(anz!=0)
      U[trans27((Ort1D)(2*i),(Ort1D)(2*j),MOrt)] = 
        U[trans27((Ort1D)(2*i),(Ort1D)(2*j),MOrt)] / anz;
    else U[trans27((Ort1D)(2*i),(Ort1D)(2*j),MOrt)] = 0.0;
  }
}

void Restrict_local_stiffness_matrix(double*  Stencil_coarse,
                                     double** Stencils_fine,
                                     bool* cid) {

  double Ucoarse[27];
  double Vcoarse[27];
  int i,j,k,l,m;
  int trans_l, trans_m;
  
  // 1.Part: Set zero
  for(i=0;i<8;++i) for(j=0;j<8;++j) Stencil_coarse[j+8*i] = 0.0;
  
  // 2.Part: Iterate
  for(i=0;i<8;++i) { // iteration for U
    Set_coarse_function(Ucoarse,(dir_sons)i,cid);
    for(j=0;j<8;++j) { // iteration for V
      Set_coarse_function(Vcoarse,(dir_sons)j,cid);
      for(k=0;k<8;++k) { // iteration over cells
        for(l=0;l<8;++l) { // iteration over corners of cell for U
          trans_l = trans_cellP_27((dir_sons)k,(dir_sons)l);
          for(m=0;m<8;++m) { // iteration over corners of cell for V
            trans_m = trans_cellP_27((dir_sons)k,(dir_sons)m);
            Stencil_coarse[j+8*i] = Stencil_coarse[j+8*i] +
              Ucoarse[trans_l] * Vcoarse[trans_m] * 
              Stencils_fine[k][m+8*l];
          }
        }
      }
    }
  }

}

---