src/expde/domain/sca_dom.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
//

// ------------------------------------------------------------
//
// sca_dom.cc
//
// ------------------------------------------------------------


// Include level 0:
#ifdef COMP_GNUOLD
 #include <iostream.h>
 #include <fstream.h>
 #include <strstream.h>
 #include <math.h>
 #define Problemzeile istrstream line(zeile,100);
#else
 #include <iostream>
 #include <fstream>
 #include <sstream> 
 #include <cmath>
 #define Problemzeile std::istringstream line(zeile,std::ios::in);
#endif 





// 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/sca_dom.h"

Scat_domain::Scat_domain(int number) {
  sample = new D3vector[number];
  normal = new D3vector[number];
  punkte = number;
}

Six_corner_pyramid::Six_corner_pyramid(
                    D3vector top1, D3vector top2, D3vector top3, 
                    D3vector top4, D3vector top5, D3vector top6,
                    D3vector down1, D3vector down2, D3vector down3, 
                    D3vector down4, D3vector down5, D3vector down6) :
  Scat_domain(8) {
  double normV;

  // 1.Step Calc border
  bordermax = top1;
  bordermax = MAX(top1,bordermax);
  bordermax = MAX(top2,bordermax);
  bordermax = MAX(top3,bordermax);
  bordermax = MAX(top4,bordermax);
  bordermax = MAX(top5,bordermax);
  bordermax = MAX(top6,bordermax);

  bordermax = MAX(down1,bordermax);
  bordermax = MAX(down2,bordermax);
  bordermax = MAX(down3,bordermax);
  bordermax = MAX(down4,bordermax);
  bordermax = MAX(down5,bordermax);
  bordermax = MAX(down6,bordermax);

  bordermin = top1;
  bordermin = MIN(top1,bordermin);
  bordermin = MIN(top2,bordermin);
  bordermin = MIN(top3,bordermin);
  bordermin = MIN(top4,bordermin);
  bordermin = MIN(top5,bordermin);
  bordermin = MIN(top6,bordermin);

  bordermin = MIN(down1,bordermin);
  bordermin = MIN(down2,bordermin);
  bordermin = MIN(down3,bordermin);
  bordermin = MIN(down4,bordermin);
  bordermin = MIN(down5,bordermin);
  bordermin = MIN(down6,bordermin);


  // 2.Step Calc bounding cube
  H    = L_infty(bordermax-bordermin)*Bounding_cube_enlargement;
  VecH = (bordermax-bordermin)*Bounding_cube_enlargement;
  A = (bordermin+bordermax - VecH)/2;
  
  // 3.Step Calc normal and sample
  // top
  sample[0] = top1;
  normal[0] = cross_product(top2-top1,   top4-top1);

  // down
  sample[1] = down1;
  normal[1] = cross_product(down4-down1, down2-down1);

  // sides
  sample[2] = top1;
  normal[2] = cross_product(down2-top1,   top2-top1);

  sample[3] = top2;
  normal[3] = cross_product(down3-top2,   top3-top2);

  sample[4] = top3;
  normal[4] = cross_product(down4-top3,   top4-top3);

  sample[5] = top4;
  normal[5] = cross_product(down5-top4,   top5-top4);

  sample[6] = top5;
  normal[6] = cross_product(down6-top5,   top6-top5);

  sample[7] = top6;
  normal[7] = cross_product(down1-top6,   top1-top6);


  for(int i=0;i<8;++i) {
    normV = norm(normal[i]);
    normal[i] = normal[i] / normV;
  }

}



Scat_domain::Scat_domain(ifstream& Datei)
{
    // 1.Step: Count number of samples
    char zeile[100];
    bordermax = D3vector(-9.999e9);
    bordermin = D3vector(9.999e9);
  
    punkte = 0;
    do {
      Datei.getline(zeile,100);
      if ((zeile[0]==' ' || zeile[0]=='-') &&
          (zeile[1]>='0' && zeile[1]<='9')) punkte++;
    }
    while ( !Datei.eof() );
    // 2.Step: initialze storage
    sample = new D3vector[punkte];
    normal = new D3vector[punkte];
    // 3.Step: read samples
    Datei.clear();  // clear EOF-bit
    Datei.seekg(0); // reset readpointer to beginning
    punkte = 0;
    char c; // dummy

    // ada std::_Ios_Openmode openMode = std::ios_base::in;;

    do {
      Datei.getline( zeile,100 );
      if ((zeile[0]==' ' || zeile[0]=='-') && (zeile[1]>='0' && zeile[1]<='9'))
        {
          Problemzeile
          line >> sample[punkte].x >> c >> sample[punkte].y >> c 
               >> sample[punkte].z >> c;
          line >> normal[punkte].x >> c >> normal[punkte].y >> c 
               >> normal[punkte].z;
          bordermax.x = 
            (bordermax.x<sample[punkte].x)?sample[punkte].x:bordermax.x;
          bordermax.y = 
            (bordermax.y<sample[punkte].y)?sample[punkte].y:bordermax.y;
          bordermax.z = 
            (bordermax.z<sample[punkte].z)?sample[punkte].z:bordermax.z;
          bordermin.x =
            (bordermin.x>sample[punkte].x)?sample[punkte].x:bordermin.x;
          bordermin.y = 
            (bordermin.y>sample[punkte].y)?sample[punkte].y:bordermin.y;
          bordermin.z = 
            (bordermin.z>sample[punkte].z)?sample[punkte].z:bordermin.z;
          punkte++;
        }
    }
    while ( !Datei.eof() );
    // 4.Step: Calc bounding cube
    H = L_infty(bordermax-bordermin)*Bounding_cube_enlargement;
    A = (bordermin+bordermax-D3vector(H))/2;
    
    bordermin = bordermin -
      D3vector(H * (1.0 - 1.0/Bounding_cube_enlargement));
    bordermax = bordermax +
      D3vector(H * (1.0 - 1.0/Bounding_cube_enlargement));

    // 5.Step: Set n_uniform
    n_uniform = 3;  // ????
    return;
};

bool Scat_domain::point_in_domain(D3vector V)
{
  if(!(bordermin<V && V<bordermax)) return false;
    D3vector delta;
    double dist;

    /* version zahn
    delta=V-sample[0];
    d2min=product(delta,delta);
    dist=product(delta,normal[0]);
    for (int i=1;i<punkte;i++) {
      delta=V-sample[i];
      d2=product(delta,delta);
      if (d2<d2min) {
        d2min=d2;
        dist=product(delta,normal[i]);
      }
    }
    return (dist<0.0);
    */

    for(int i=0;i<punkte;i++) {
      delta=V-sample[i];
      dist=product(delta,normal[i]);
      if(dist>=0.0) return false;
    }
    return true;
}

double Scat_domain::distance(D3vector V, dir_3D d) {
  if(!(bordermin<V && V<bordermax)) return -100.0;

  D3vector dir_vector,delta,normal_vec, normal_vec2;
  switch (d)
    {
    case Wdir: dir_vector=D3vector(-1,0,0); break;
    case Edir: dir_vector=D3vector(1,0,0); break;
    case Sdir: dir_vector=D3vector(0,-1,0); break;
    case Ndir: dir_vector=D3vector(0,1,0); break;
    case Ddir: dir_vector=D3vector(0,0,-1); break;
    default: dir_vector=D3vector(0,0,1);
    }

  double dis_R, dis_R_new, dist, prod_dir_nor;
  dis_R = 4.0*H;
  for(int i=0;i<punkte;i++) {
    prod_dir_nor = product(dir_vector,normal[i]);
    if(prod_dir_nor > 0) {
      delta=V-sample[i];
      dist=product(delta,normal[i]);
      dis_R_new = -dist/prod_dir_nor;

      if(dis_R_new<dis_R)
        dis_R = dis_R_new;
    }
  }
  return dis_R;
  

  /* Version zahn
  double dist, d2, d2min, d2min_in_dom, dist_in_dom;

  dist = -100.0;
  dist_in_dom = 10;
  
  delta=V-sample[0];
  d2min=product(delta,delta);
  
  d2min = 10e10;
  d2min_in_dom = 10e10;
  for (int i=0;i<punkte;i++) {
    delta=V-sample[i];
    d2=product(delta,delta);
    normal_vec2=normal[i];
    if(d2<d2min_in_dom) {
      d2min_in_dom=d2;
      dist_in_dom=product(delta,normal[i]);
    }
    //      if(d2<d2min && product(dir_vector,normal_vec2)>0.0) {
    //      if(d2<d2min && product(dir_vector,normal_vec2)>0.1) {
    if(d2<d2min && product(dir_vector,normal_vec2)>0.01) {
      d2min=d2;
      normal_vec = normal_vec2;
      dist=product(delta,normal[i]);
    }
  }
  
  //    return 0.0;
  
  if(developer_version) {
    if(product(dir_vector,normal_vec)<0.0
       && point_in_domain(V)) {
      cout << " \n Mistake wrong direction" << dist << endl;
    }
  }
  
  if(dist_in_dom < 0.0)
    return -dist/product(dir_vector,normal_vec);
  else return -100.0;
  */
}

---