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

// ------------------------------------------------------------
//
// domain.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 2:
#include "../basic/basic.h"

// Include level 3:
#include "domain.h"

#include "d_exam.h"


//  SQUARE:      [0,1]^3
// ----------------------------------------------

PointtypeD All_interior(D3vector V) {
  if(V.x < 1.0 && V.x > 0.0 &&
     V.y < 1.0 && V.y > 0.0 &&
     V.z < 1.0 && V.z > 0.0)
    return interiorD;
  else return exteriorD; 
};


double dis_Square(D3vector V, dir_3D d) {
  if(V.x < 1.0 && V.x > 0.0 &&
     V.y < 1.0 && V.y > 0.0 &&
     V.z < 1.0 && V.z > 0.0) {
    if(d==Tdir) return 1.0 - V.z;
    if(d==Ddir) return V.z;
    
    if(d==Edir) return 1.0 - V.x;
    if(d==Wdir) return V.x;
    
    if(d==Ndir) return 1.0 - V.y;
    if(d==Sdir) return V.y;
  }
  return 0.0;
}


//  BALL:         Radius: 0.5,  center: (0.5,0.5,0.5)
// ---------------------------------------------------

#define radius_ball 0.5

PointtypeD Poi_Ball(D3vector V) { 
  V.x = V.x - 0.5;
  V.y = V.y - 0.5;
  V.z = V.z - 0.5;
  
    if(sqrt(V.x*V.x+V.y*V.y+V.z*V.z)<=radius_ball) return interiorD;
    else return exteriorD;
}

double dis_Ball(D3vector V, dir_3D d) {
  static double x,y,z; 
  x = V.x - 0.5;
  y = V.y - 0.5;
  z = V.z - 0.5;
  if(sqrt(x*x+y*y+z*z) < radius_ball) {
    if(d==Tdir) return sqrt(radius_ball*radius_ball-x*x-y*y) - z;
    if(d==Ddir) return sqrt(radius_ball*radius_ball-x*x-y*y) + z;

    if(d==Edir) return sqrt(radius_ball*radius_ball-y*y-z*z) - x;
    if(d==Wdir) return sqrt(radius_ball*radius_ball-y*y-z*z) + x;

    if(d==Ndir) return sqrt(radius_ball*radius_ball-x*x-z*z) - y;
    if(d==Sdir) return sqrt(radius_ball*radius_ball-x*x-z*z) + y;
  }
  return 0.0;
}


//  Cylinder:         Radius: 0.5,  center: (0.5,0.5,0.5)
// ---------------------------------------------------

#define radius_cyl 0.5

PointtypeD Poi_Cylinder(D3vector V) { 
  V.x = V.x - 0.5;
  V.y = V.y - 0.5;
  
  if(sqrt(V.x*V.x+V.y*V.y)<=radius_cyl
     && V.z>=0.0 && V.z<=1.0) return interiorD;
  else return exteriorD;
}

double dis_Cylinder(D3vector V, dir_3D d) {
  static double x,y; 
  x = V.x - 0.5;
  y = V.y - 0.5;
  if(sqrt(x*x+y*y)<=radius_cyl
     && V.z>=0.0 && V.z<=1.0) {
    if(d==Tdir) return 1.0 - V.z;
    if(d==Ddir) return V.z;
    
    if(d==Edir) return sqrt(radius_cyl*radius_cyl-y*y) - x;
    if(d==Wdir) return sqrt(radius_cyl*radius_cyl-y*y) + x;

    if(d==Ndir) return sqrt(radius_cyl*radius_cyl-x*x) - y;
    if(d==Sdir) return sqrt(radius_cyl*radius_cyl-x*x) + y;
  }
  return 0.0;
}



//  Skew_cylinder:  radius: r_down, r_top,   center: (0.0,0.0,0.5)
// ------------------------------------------------------------
//  R = R_down * V.z + R_top * (1-V.z)

PointtypeD Poi_Skew_cylinder(D3vector V, double R_down, double R_top) { 
  double radius;
  radius = R_top * V.z + R_down * (1.0-V.z);
  
  if(sqrt(V.x*V.x+V.y*V.y)<=radius
     && V.z>=0.0 && V.z<=1.0) return interiorD;
  else return exteriorD;
}

double dis_Skew_cylinder(D3vector V, dir_3D d, double R_down, double R_top) {
  double radius, r_point;
  radius = R_top * V.z + R_down * (1.0-V.z);
  r_point = sqrt(V.x*V.x+V.y*V.y);
  if(r_point<=radius
     && V.z>=0.0 && V.z<=1.0) {
    if(d==Ddir) {
      if(r_point<R_down) {
        return V.z;
      }
      else {
        return (r_point + R_top) / (R_down - R_top);
      }
    }
    if(d==Tdir) {
      if(r_point<R_top) {
        return 1.0-V.z;
      }
      else {
        return 1.0-(r_point + R_top) / (R_down - R_top);
      }
    }

    if(d==Edir) return sqrt(radius*radius-V.y*V.y) - V.x;
    if(d==Wdir) return sqrt(radius*radius-V.y*V.y) + V.x;

    if(d==Ndir) return sqrt(radius*radius-V.x*V.x) - V.y;
    if(d==Sdir) return sqrt(radius*radius-V.x*V.x) + V.y;
  }
  return 0.0;
}

//  Double_cylinder:  radius: R_down = r_in, 
//                            R_top  = r_out,   center: (0.0,0.0,0.5)
// ------------------------------------------------------------
//  R = R_down * V.z + R_top * (1-V.z)

PointtypeD Poi_Double_cylinder(D3vector V, double r_in, double r_out) { 
  double radius;
  radius = sqrt(V.x*V.x+V.y*V.y);
  
  if(V.z>=0.0 && V.z<=1.0) {
    if(r_in<=radius && radius<=r_out)
      return interiorD;
  }
  return exteriorD;
}

double dis_Double_cylinder(D3vector V, dir_3D d, double r_in, double r_out) {
  double radius;
  radius = sqrt(V.x*V.x+V.y*V.y);
  if(r_in<=radius && radius<=r_out
     && V.z>=0.0 && V.z<=1.0) {
    if(d==Ddir) {
      return V.z;
    }
    if(d==Tdir) {
      return 1.0-V.z;
    }

    if(d==Edir) {
      if(V.y < -r_in || V.y > r_in || V.x > 0.0)
        return sqrt(r_out*r_out-V.y*V.y) - V.x;
      else
        return - sqrt(r_in*r_in-V.y*V.y) - V.x;
    }
    if(d==Wdir) {
      if(V.y < -r_in || V.y > r_in || V.x < 0.0)
        return sqrt(r_out*r_out-V.y*V.y) + V.x;
      else
        return - sqrt(r_in*r_in-V.y*V.y) + V.x;
    }
    if(d==Ndir) {
      if(V.x < -r_in || V.x > r_in || V.y > 0.0)
        return sqrt(r_out*r_out-V.x*V.x) - V.y;
      else
        return - sqrt(r_in*r_in-V.x*V.x) - V.y;
    }
    if(d==Sdir) {
      if(V.x < -r_in || V.x > r_in || V.y < 0.0)
        return sqrt(r_out*r_out-V.x*V.x) + V.y;
      else
        return - sqrt(r_in*r_in-V.x*V.x) + V.y;
    }
  }
  return 0.0;
}



//  SCREW SQUARE:  angles:               alpha_up, alpha_down
// -------------------------------------------------------------


PointtypeD Poi_Skew_square(D3vector V, double R_down, double R_top) { 
  if(0<=V.x     && V.x<=1.0   &&
     0<=V.y     && V.y<=1.0   &&
     V.z>=my_tan(R_down)*V.x    &&
     V.z<=1.0+my_tan(R_top)*V.x) return interiorD;
  else return exteriorD;
}

double dis_Skew_square(D3vector V, dir_3D d, double R_down, double R_top) {
  if(0<=V.x     && V.x<=1.0   &&
     0<=V.y     && V.y<=1.0   &&
     V.z>=my_tan(R_down)*V.x    &&
     V.z<=1.0+my_tan(R_top)*V.x) {
    if(d==Tdir) return 1.0+my_tan(R_top)*V.x - V.z;
    if(d==Ddir) return V.z-my_tan(R_down)*V.x;
    
    if(d==Edir) {
      if(V.z>1.0+my_tan(R_top)) {
        return (V.z-1.0)/my_tan(R_top)-V.x;
      }
      if(V.z<my_tan(R_down)) {
        return V.z/my_tan(R_down)-V.x;
      }
      return 1.0 - V.x;
    }
    if(d==Wdir) {
      if(V.z>1.0) {
        return V.x-(V.z-1.0)/my_tan(R_top);
      }
      if(V.z<0.0) {
        return V.x-V.z/my_tan(R_down);
      }
      return V.x;
    }
    if(d==Ndir) return 1.0 - V.y;
    if(d==Sdir) return V.y;
  }
  return 0.0;
}



//  SCREW SQUARE: DTLR:  angles:            alpha_up,   alpha_down
//                                         alpha_left, alpha_right
// ----------------------------------------------------------------

/*
PointtypeD Poi_Skew_squareDTLR(D3vector V, 
                               double R_down, double R_top,
                               double R_left, double R_right) { 
  if(0<=V.x     && V.x<=1.0     &&
     V.y>=my_tan(R_left)*V.x       &&
     V.y<=1.0+my_tan(R_right)*V.x  &&
     V.z>=my_tan(R_down)*V.x       &&
     V.z<=1.0+my_tan(R_top)*V.x) return interiorD;
  else return exteriorD;
}

double dis_Skew_squareDTLR(D3vector V, dir_3D d, 
                           double R_down, double R_top,
                           double R_left, double R_right) {
  double disA;
  double disB;

  if(0<=V.x     && V.x<=1.0     &&
     V.y>=my_tan(R_left)*V.x       &&
     V.y<=1.0+my_tan(R_right)*V.x &&
     V.z>=my_tan(R_down)*V.x       &&
     V.z<=1.0+my_tan(R_top)*V.x) {
    if(d==Tdir) return 1.0+my_tan(R_top)*V.x - V.z;
    if(d==Ddir) return V.z-my_tan(R_down)*V.x;
    
    if(d==Edir) {
      // for z-x
      if(V.z>1.0+my_tan(R_top)) {
        disA = (V.z-1.0)/my_tan(R_top)-V.x;
      }
      else {
        if(V.z<my_tan(R_down)) {
          disA = V.z/my_tan(R_down)-V.x;
        }
        else 
          disA = 1.0 - V.x;
      }
      // for y-x
      if(V.y>1.0+my_tan(R_right)) {
        disB = (V.y-1.0)/my_tan(R_right)-V.x;
      }
      else {
        if(V.y<my_tan(R_left)) {
          disB = V.y/my_tan(R_left)-V.x;
        }
        else 
          disB = 1.0 - V.x;
      }
      //      return disB;
      return MIN(disA,disB);
    }
    if(d==Wdir) {
      // for z-x
      if(V.z>1.0) {
        disA =  V.x-(V.z-1.0)/my_tan(R_top);
      } 
      else {
        if(V.z<0.0) {
          disA =  V.x-V.z/my_tan(R_down);
        }
        else disA =  V.x;
      }
      // for y-x
      if(V.y>1.0) {
        disB =  V.x-(V.y-1.0)/my_tan(R_right);
      }
      else {
        if(V.y<0.0) {
          disB =  V.x-V.y/my_tan(R_left);
        }
        else disB =  V.x;
      }
      return MIN(disA,disB);
    }
    if(d==Ndir) return 1.0+my_tan(R_right)*V.x - V.y;
    if(d==Sdir) return V.y-my_tan(R_left)*V.x;
  }
  return 0.0;
}
*/

PointtypeD Poi_Skew_squareDTLR(D3vector V, 
                               double R_down, double R_top,
                               double R_left, double R_right) { 
  if(0<=V.z     && V.z<=1.0     &&
     V.y>=my_tan(R_left)*V.z       &&
     V.y<=1.0+my_tan(R_right)*V.z  &&
     V.x>=my_tan(R_down)*V.z       &&
     V.x<=1.0+my_tan(R_top)*V.z) return interiorD;
  else return exteriorD;
}

double dis_Skew_squareDTLR(D3vector V, dir_3D d, 
                           double R_down, double R_top,
                           double R_left, double R_right) {
  double disA;
  double disB;

  if(0<=V.z     && V.z<=1.0     &&
     V.y>=my_tan(R_left)*V.z       &&
     V.y<=1.0+my_tan(R_right)*V.z &&
     V.x>=my_tan(R_down)*V.z       &&
     V.x<=1.0+my_tan(R_top)*V.z) {
    if(d==Edir) return 1.0+my_tan(R_top)*V.z - V.x;
    if(d==Wdir) return V.x-my_tan(R_down)*V.z;
    
    if(d==Tdir) {
      // for z-x
      if(V.x>1.0+my_tan(R_top)) {
        disA = (V.x-1.0)/my_tan(R_top)-V.z;
      }
      else {
        if(V.x<my_tan(R_down)) {
          disA = V.x/my_tan(R_down)-V.z;
        }
        else 
          disA = 1.0 - V.z;
      }
      // for y-x
      if(V.y>1.0+my_tan(R_right)) {
        disB = (V.y-1.0)/my_tan(R_right)-V.z;
      }
      else {
        if(V.y<my_tan(R_left)) {
          disB = V.y/my_tan(R_left)-V.z;
        }
        else 
          disB = 1.0 - V.z;
      }
      //      return disB;
      return MIN(disA,disB);
    }
    if(d==Ddir) {
      // for z-x
      if(V.x>1.0) {
        disA =  V.z-(V.x-1.0)/my_tan(R_top);
      } 
      else {
        if(V.x<0.0) {
          disA =  V.z-V.x/my_tan(R_down);
        }
        else disA =  V.z;
      }
      // for y-x
      if(V.y>1.0) {
        disB =  V.z-(V.y-1.0)/my_tan(R_right);
      }
      else {
        if(V.y<0.0) {
          disB =  V.z-V.y/my_tan(R_left);
        }
        else disB =  V.z;
      }
      return MIN(disA,disB);
    }
    if(d==Ndir) return 1.0+my_tan(R_right)*V.z - V.y;
    if(d==Sdir) return V.y-my_tan(R_left)*V.z;
  }
  return 0.0;
}




//  periodic Cylinder:         Radius: 0.5,  center: (0.5,0.5,0.5)
// ---------------------------------------------------


PointtypeD Poi_periodic_Cylinder(D3vector V) { 
  V.x = V.x - 0.5;
  V.y = V.y - 0.5;
  
  if(sqrt(V.x*V.x+V.y*V.y)<=radius_cyl) return interiorD;
  else return exteriorD;
}

double dis_periodic_Cylinder(D3vector V, dir_3D d) {
  static double x,y; 
  x = V.x - 0.5;
  y = V.y - 0.5;
  if(sqrt(x*x+y*y)<=radius_cyl
     && V.z>=0.0 && V.z<=1.0) {
    if(d==Tdir || d==Ddir) {
      cout << " error in periodic cylinder " << endl;
      return 10.0e10;
    } 
   
    if(d==Edir) return sqrt(radius_cyl*radius_cyl-y*y) - x;
    if(d==Wdir) return sqrt(radius_cyl*radius_cyl-y*y) + x;

    if(d==Ndir) return sqrt(radius_cyl*radius_cyl-x*x) - y;
    if(d==Sdir) return sqrt(radius_cyl*radius_cyl-x*x) + y;
  }
  return 0.0;
}

---