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

// ------------------------------------------------------------
//
// input.cc
//
// ------------------------------------------------------------



// Include level 0:
#ifdef COMP_GNUOLD
 #include <iostream.h>
 #include <fstream.h>
 #include <math.h>
 #include <time.h>
#else
 #include <iostream>
 #include <fstream>
 #include <cmath>
 #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 "gpar.h"
#include "parallel.h"
#include "mgcoeff.h"
#include "sto_man.h"
#include "gridbase.h"
#include "grid.h"
#include "input.h"



Input_data_object::Input_data_object (D3vector* coordinates,
                                      double* value,
                                      int num_points,
                                      double mesh_x,
                                      double mesh_y,
                                      double mesh_z) {
  int i;
  
  // 1. set 
  ort   = coordinates;
  werte = value;
  anz_punkte = num_points;
  h_x = mesh_x;
  h_y = mesh_y;
  h_z = mesh_z;

  epsx = h_x * factor_eps_input;
  epsy = h_y * factor_eps_input;
  epsz = h_z * factor_eps_input;

  // 2. hashtable
  hashtable_input_leng  = Find_next_prime(2*num_points);
  hashtable_input_start = 
    new Point_hashtable_input*[hashtable_input_leng];
  for(i=0;i<hashtable_input_leng;++i)
    hashtable_input_start[i] = NULL;
  hashtable_input_occ=0;

  // 3. find minimum
  org_x = coordinates[0].x;
  org_y = coordinates[0].y;
  org_z = coordinates[0].z;
  for(i=0;i<num_points;++i) {
    org_x = MIN(coordinates[0].x,org_x);
    org_y = MIN(coordinates[0].y,org_y);
    org_z = MIN(coordinates[0].z,org_z);
  }

  // 4. add points
  for(i=0;i<num_points;++i)
    Add_point(coordinates[i].x,coordinates[i].y,coordinates[i].z,value[i]);

}

void Input_data_object::Add_point(double x, double y, double z,
                                  double value) {
  Point_hashtable_input* point;
  int Nx, Ny, Nz;

  Nx = (int)((x+h_x/4.0-org_x)/h_x);
  Ny = (int)((y+h_y/4.0-org_y)/h_y);
  Nz = (int)((z+h_z/4.0-org_z)/h_z);
  

  point = hashtable_input_start
    [hashtable_input_function(Nx,Ny,Nz,hashtable_input_leng)];
  if(point==NULL) {
    hashtable_input_occ++;
    point = new Point_hashtable_input(Nx,Ny,Nz,value);
    hashtable_input_start[hashtable_input_function(Nx,Ny,Nz,
                                                   hashtable_input_leng)] = 
      point;
  }
  else {
    while(point->next!=NULL && 
          (point->nx!=Nx || point->ny!=Ny || point->nz!=Nz))
      point = point->next;
    if(point->nx!=Nx || point->ny!=Ny || point->nz!=Nz) {
      hashtable_input_occ++;
      point->next = new Point_hashtable_input(Nx,Ny,Nz,value);
    }
  }
}


double Input_data_object::Give_value(D3vector coordinates) {
  int Nx, Ny, Nz;
  Point_hashtable_input* point;

  Nx = (int)((coordinates.x+epsx-org_x)/h_x);
  Ny = (int)((coordinates.y+epsy-org_y)/h_y);
  Nz = (int)((coordinates.z+epsz-org_z)/h_z);
  
  point = 
    hashtable_input_start[hashtable_input_function(Nx,Ny,Nz,
                                     hashtable_input_leng)];
  while(point!=NULL) {
    if(point->Is_point(Nx,Ny,Nz)) return point->wert;
    point = point->next;
  }
  return 0.0;
}

bool Input_data_object::Exists_point(D3vector coordinates) {
  int Nx, Ny, Nz;
  Point_hashtable_input* point;

  Nx = (int)((coordinates.x+epsx-org_x)/h_x);
  Ny = (int)((coordinates.y+epsy-org_y)/h_y);
  Nz = (int)((coordinates.z+epsz-org_z)/h_z);
  
  point = 
    hashtable_input_start[hashtable_input_function(Nx,Ny,Nz,
                                     hashtable_input_leng)];
  while(point!=NULL) {
    if(point->Is_point(Nx,Ny,Nz)) return true;
    point = point->next;
  }
  return false;
}

D3vector Input_data_object::Staggered_coordinate(D3vector coordinates) {
  return D3vector(((int)((coordinates.x+epsx-org_x)/h_x))*h_x+org_x,
                  ((int)((coordinates.y+epsy-org_y)/h_y))*h_y+org_y,
                  ((int)((coordinates.z+epsz-org_z)/h_z))*h_z+org_z);
}


double Input_data_object::Interpolate_value(D3vector coordinates) {
  D3vector lambda;
  int i,j,k,l;
  double sum;
  int anz;

  if(Exists_point(coordinates) &&
     Exists_point(coordinates+D3vector(h_x,0.0,0.0)) &&
     Exists_point(coordinates+D3vector(0.0,h_y,0.0)) &&
     Exists_point(coordinates+D3vector(0.0,0.0,h_z)) &&
     Exists_point(coordinates+D3vector(0.0,h_y,h_z)) &&
     Exists_point(coordinates+D3vector(h_x,0.0,h_z)) &&
     Exists_point(coordinates+D3vector(h_x,h_y,0.0)) &&
     Exists_point(coordinates+D3vector(h_x,h_y,h_z))) {
    lambda = (coordinates - Staggered_coordinate(coordinates))/
      D3vector(h_x,h_y,h_z);
    
    sum = 0.0;
    for(i=0;i<2;++i) for(j=0;j<2;++j) for(k=0;k<2;++k) {
      sum = sum +
        Give_value(coordinates+D3vector(i*h_x,j*h_y,k*h_z))
        * ((1.0*(1-i) - lambda.x * (1-2*i)))
        * ((1.0*(1-j) - lambda.y * (1-2*j)))
        * ((1.0*(1-k) - lambda.z * (1-2*k)));
    }
    return sum;
  }
  else {
    for(l=1;l<search_length;++l) {
      anz = 0;
      sum = 0.0;
      for(i=0;i<2;++i) for(j=0;j<2;++j) for(k=0;k<2;++k) {
        if(Exists_point(coordinates+D3vector(h_x*((1-l)+i*(2*l-1)),
                                             h_y*((1-l)+j*(2*l-1)),
                                             h_z*((1-l)+k*(2*l-1))))) {
          sum = sum +
            Give_value(coordinates+D3vector(h_x*((1-l)+i*(2*l-1)),
                                            h_y*((1-l)+j*(2*l-1)),
                                            h_z*((1-l)+k*(2*l-1))));
          anz =  anz + 1;
        }
      }
      if(anz>0) return sum / (double)anz;
    }
  }
  return 0.0;
}

     

---