src/expde/extemp/partint.cc

Go to the documentation of this file.

#include "partint.h"


bool In_tet(dir_sons Tet[4], D3vector corn[8], D3vector V, double h) {

  const double EPS = -1.0e-5*h*h*h*h;

  if(Calc_bary_weights(corn[Tet[0]],corn[Tet[1]],corn[Tet[2]],corn[Tet[3]],V,0)<=EPS)
    return false;
  if(Calc_bary_weights(corn[Tet[0]],corn[Tet[1]],corn[Tet[2]],corn[Tet[3]],V,1)<=EPS)
    return false;
  if(Calc_bary_weights(corn[Tet[0]],corn[Tet[1]],corn[Tet[2]],corn[Tet[3]],V,2)<=EPS)
    return false;
  if(Calc_bary_weights(corn[Tet[0]],corn[Tet[1]],corn[Tet[2]],corn[Tet[3]],V,3)<=EPS)
    return false;
  return true;
}

bool In_tetbo(Tetraeder_storage* tet, D3vector corn[28], D3vector V) {
  if(Calc_bary_weights(corn[tet->N0()],corn[tet->N1()],
                       corn[tet->N2()],corn[tet->N3()],
                       V,0)<0.0)
    return false;
  if(Calc_bary_weights(corn[tet->N0()],corn[tet->N1()],
                       corn[tet->N2()],corn[tet->N3()],
                       V,1)<0.0)
    return false;
  if(Calc_bary_weights(corn[tet->N0()],corn[tet->N1()],
                       corn[tet->N2()],corn[tet->N3()],
                       V,2)<0.0)
    return false;
  if(Calc_bary_weights(corn[tet->N0()],corn[tet->N1()],
                       corn[tet->N2()],corn[tet->N3()],
                       V,3)<0.0)
    return false;
  return true;
}

void Set_tets(dir_sons Tets[6][4]) {
  // 0. Tet: WNT, WND, WST, EST 
  Tets[0][0] = WNTd;
  Tets[0][1] = WNDd;
  Tets[0][2] = WSTd; 
  Tets[0][3] = ESTd;     

  // 1. Tet: EST, WND, WST, ESD     
  Tets[1][0] = ESTd;
  Tets[1][1] = WNDd;
  Tets[1][2] = WSTd; 
  Tets[1][3] = ESDd;     

  // 2. Tet: WND, WSD, WST, ESD        
  Tets[2][0] = WNDd;
  Tets[2][1] = WSDd;
  Tets[2][2] = WSTd; 
  Tets[2][3] = ESDd;     

  // 3. Tet: EST, WND, ESD, END
  Tets[3][0] = ESTd;
  Tets[3][1] = WNDd;
  Tets[3][2] = ESDd; 
  Tets[3][3] = ENDd;     

  // 4. Tet: ENT, WNT, EST, END
  Tets[4][0] = ENTd;
  Tets[4][1] = WNTd;
  Tets[4][2] = ESTd; 
  Tets[4][3] = ENDd;     

  // 5. Tet: WNT, WND, EST, END
  Tets[5][0] = WNTd;
  Tets[5][1] = WNDd;
  Tets[5][2] = ESTd; 
  Tets[5][3] = ENDd;     
}



void partHelm_FE(Variable *f, double x, double y, double z, double q) {
  Grid* grid;
  grid = f->Give_grid();
  Index3D I,J;
  D3vector V;
  const int max_lev = grid->Max_level();
  double h;
  dir_sons Tets[6][4];
  int num;
  
  h   = grid->Give_finest_mesh_size();
  V = grid->back_transform_coord(D3vector(x,y,z));
  I = back_coordinate(V,max_lev);
  Set_tets(Tets);

  D3vector corn[28];
  bool found;
  Tetraeder_storage* found_tet;

  for(unsigned int i=0;i<8;++i) {
    corn[i] = D3vector(x1DCoord((dir_sons)i),
                       y1DCoord((dir_sons)i),
                       z1DCoord((dir_sons)i));
  }
  
  if(grid->Give_cell_typ(I)==int_cell) {
    /*
    // Gehe zu 8 Ecken der Zelle 
    for(unsigned int i=0;i<8;++i) {
      J=I.neighbour((dir_sons)i);
      
      grid->Give_variable(J,max_lev)[f->Number_variable()] += 1.0/8.0 * q;
    }
    */
    V = (D3vector(x,y,z) - grid->transform_coord(I.neighbour(WSDd))) / h;

    found=false;
    int ii=0;
    for(;ii<6&&(found==false);++ii) {
      found = In_tet(Tets[ii],corn,V,h);
    }
    if(found==false) cout << " error in partHelm_FE found tet! " << endl;
    --ii;

    for(unsigned int i=0;i<4;++i) {
      J=I.neighbour(Tets[ii][i]);
      
      grid->Give_variable(J,max_lev)[f->Number_variable()] +=
        q * Calc_bary_weights(corn[Tets[ii][0]],corn[Tets[ii][1]],
                              corn[Tets[ii][2]],corn[Tets[ii][3]],V,i);
    }
  } else if(grid->Give_cell_typ(I)==fine_bo_cell) {
    /*
    int count = 0;    
    for(unsigned int i=0;i<8;++i) {
      J=I.neighbour((dir_sons)i);
      if (grid->Give_type(J) == nearb)
        count++;
    }
    for(unsigned int i=0;i<8;++i) {
      J=I.neighbour((dir_sons)i);
      if (grid->Give_type(J) == nearb) {
        grid->Give_variable(J,max_lev)[f->Number_variable()] += 1.0/(double)count * q;
      }
    }
    */
    
    BoCell* boc = grid->Give_Bo_cell(I);
    for(num=0;num<boc->Give_total_number_points();++num) {
      corn[num] = boc->coord(num , boc->hhh) / h;
    }

    found_tet=NULL;    
    Tetraeder_storage* tet;
    for(tet=boc->Give_tets();tet!=NULL && found_tet==NULL;tet=tet->Next()) {
      if(In_tetbo(tet,corn,V)) {
        found_tet = tet;
      }
    }
    if(found_tet==NULL) {
      int count = 0;    
      for(unsigned int i=0;i<8;++i) {
        J=I.neighbour((dir_sons)i);
        if (grid->Give_type(J) == nearb)
          count++;
      }
      for(unsigned int i=0;i<8;++i) {
        J=I.neighbour((dir_sons)i);
        if (grid->Give_type(J) == nearb) {
          grid->Give_variable(J,max_lev)[f->Number_variable()] += 1.0/(double)count * q;
        }
      }
    }
    else {      
      for(unsigned int i=0;i<4;++i) {
          if(i==0) num = found_tet->N0();
          if(i==1) num = found_tet->N1();
          if(i==2) num = found_tet->N2();
          if(i==3) num = found_tet->N3();
          boc->vars[num][f->Number_variable()] +=
              q * Calc_bary_weights(corn[found_tet->N0()],corn[found_tet->N1()],
                  corn[found_tet->N2()],corn[found_tet->N3()],
                  V,i) / found_tet->Det();
/*
  Changed that while I was in Berkley October 2004 in order to have the
  charge density instead the charges !
  
  boc->vars[num][f->Number_variable()] +=
  q * Calc_bary_weights(corn[found_tet->N0()],corn[found_tet->N1()],
  corn[found_tet->N2()],corn[found_tet->N3()],
  V,i) / found_tet->Det() *h*h*h;
*/
      }      
    }
  }
  else {
    /*
    cout << " x: " << x
         << " y: " << y
         << " z: " << z
         << " max_lev: " <<  max_lev
         << " q: " << q
         << endl;
    
    cout << "\n I: "; 
    I.Print();
    cout << "\n coord of I: "; 
    grid->transform_coord(I).Print();
    cout << endl;
    cout << " Typ der Zelle: " << (int)grid->Give_cell_typ(I) << endl;
    */
    // cout << "\n can not be .... " << endl;
  }  
  // fehlt: Parallel aufaddieren der Summation bei GhostPunkten

  /*
  int array_num[1];
  array_num[0] = f->Number_variable();
  grid->Sum_ghost_nodes_Variable(array_num,1,grid->Max_level());
  */
}

inline double weight(double x, dir1D d, double h) {
  if(Rd==d) return x;
            return h-x;
}

void partInterp(Variable *f, double x, double y, double z, double *q) {
  Grid* grid;
  grid = f->Give_grid();
  Index3D I,J;
  D3vector V;
  const int max_lev = grid->Max_level();
  double h, hhh;

  h   = grid->Give_finest_mesh_size();
  hhh = h*h*h;

  //  aus (x,y,z) Index der Zelle class Index3D I, getCell?
  V = grid->back_transform_coord(D3vector(x,y,z));
  I = back_coordinate(V,max_lev);
  
  *q=0.0;

  if(grid->Give_cell_typ(I)==int_cell) {
    /*
    // Gehe zu 8 Ecken der Zelle 
    for(unsigned int i=0;i<8;++i) {
      J=I.neighbour((dir_sons)i);
      *q += grid->Give_variable(J,max_lev)[f->Number_variable()] / 8.0;
    }
    */
      V = D3vector(x,y,z) - grid->transform_coord(I.neighbour(WSDd));
      for(unsigned int i=0;i<8;++i) {
          J=I.neighbour((dir_sons)i);
          *q += grid->Give_variable(J,max_lev)[f->Number_variable()] *
              weight(V.x,x1DCoord((dir_sons)i),h) *
              weight(V.y,y1DCoord((dir_sons)i),h) *
              weight(V.z,z1DCoord((dir_sons)i),h) / hhh;
      }
  } 
  else if(grid->Give_cell_typ(I)==fine_bo_cell) {
      int count = 0;
      for(unsigned int i=0;i<8;++i) {
          J=I.neighbour((dir_sons)i);
          if (grid->Give_type(J) == nearb)
              count++;
      }
      for(unsigned int i=0;i<8;++i) {
          J=I.neighbour((dir_sons)i);
          if (grid->Give_type(J) == nearb) {
              *q += grid->Give_variable(J,max_lev)[f->Number_variable()] / (double)count;
          }
      }
  }
  else {
      //  cout << "\n can not be .... " << endl;
  }  
  // fehlt: Parallel aufaddieren der Summation bei GhostPunkten
}

---