src/expde/domain/domain.h

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

#ifndef DOMAIN_H_
#define DOMAIN_H_
      
// 1.:  basic class
// 2.:  new domain
// 3.:  convex domains
// 4.1:  union domains
// 4.2:  cut domains
// 5.:  move domains
// 6.:  scale domains
// 7.:  operators
// 8.:  member functions 
// 9.:  special domain fuer Andreas  
// examples: d_examp.h

// 3D unit vector
//-------------------
inline D3vector Unit_vector(dir_3D d) {
  if(d>Ndir) return D3vector(0.0,0.0,d-4);
  if(d<Sdir) return D3vector(d,0.0,0.0);
  return D3vector(0.0,d-2,0.0);
}


// 1.:  basic class

class All_Domains {
 public:
  // constructor 
  All_Domains() {}; 
  
  // size of minimal uniform grid (in levels) (only for error messages!!)
  virtual int Give_n_uniform() =0;
  // is it possible to calculate whether edge in domain?
  virtual bool calc_edge(D3vector V1, D3vector V2) =0;
  // is edge in domain?
  virtual bool edge(D3vector V1, D3vector V2) =0;
  //  is it possible to calculate distance to the boundary?
  virtual calc_dis calc_distance(D3vector V, dir_3D d) =0;
  // distance to the boundary in direction d
  virtual double distance(D3vector V, dir_3D d) =0;
  // is point in domain? This function should be used only 
  // at the beginning of the grid generation
  virtual bool point_in_domain(D3vector V) =0;

  bool Test(double x) { return 1.0; }

  // describtion of the square, in which the domain is contained
  // [a1,a1+h]x[a2,a2+h]x[a3,a3+h]
  virtual D3vector GiveA() const =0;
  virtual double   GiveH() const =0;
  virtual D3vector GiveVecH() const =0;

  // Periodoic object??
  virtual bool Is_periodic() =0;
 protected: 
  // for  square, in which the domain is contained
  D3vector A;
  double H;
  D3vector VecH;
  // size of minimal uniform grid (in levels)
  int n_uniform;

  // Periodoic object??
  bool is_periodic;

  // ...
  PointtypeD (*function)(D3vector);   // point in domain?
  double (*func_distance)(D3vector,dir_3D);  // distance to the boundary

  PointtypeD (*function_2P)(D3vector,double,double); // point in domain?
  double (*func_distance_2P)(D3vector,dir_3D,double,double); 
                                        // distance to the boundary

  PointtypeD (*function_4P)(D3vector,double,double,double,double); 
                                        // point in domain?
  double (*func_distance_4P)(D3vector,dir_3D,double,double,double,double); 
                                        // distance to the boundary
 
  All_Domains *domain_A;
  All_Domains *domain_B;

  D3vector V_move;
  D3vector V_scale;

  // for parameters:
  double R_down;
  double R_top;

  double R_left;
  double R_right;
};

// 2.:  new domains

class Domain : public All_Domains {
 public:
  // constructor
  Domain(All_Domains* dom) {
    domain_A = dom;
  }
  Domain(Domain* dom) {
    domain_A = dom;
  }
  Domain(const Domain& dom) {
    domain_A = new Domain(dom.domain_A);
  }

  // size of minimal uniform grid (in levels)
  int  Give_n_uniform() { return domain_A->Give_n_uniform(); } ;
  // is it possible to calculate whether edge in domain?
  bool calc_edge(D3vector V1, D3vector V2) { return true; }
  // is edge in domain?
  bool edge(D3vector V1, D3vector V2) { return domain_A->edge(V1,V2); };
  //  is it possible to calculate distance to the boundary?
  calc_dis calc_distance(D3vector V, dir_3D d) {return yes; }
  // distance to the boundary in direction d
  double distance(D3vector V, dir_3D d) { return domain_A->distance(V,d); };
  // is point in domain?
  bool point_in_domain(D3vector V) { return domain_A->point_in_domain(V); }

  // Periodoic object??
  bool Is_periodic() { return domain_A->Is_periodic(); }

  // describtion of the square, in which the domain is contained
  D3vector GiveA()    const { return domain_A->GiveA(); }
  double   GiveH()    const { return domain_A->GiveH(); }
  D3vector GiveVecH() const { return domain_A->GiveVecH(); }
};

// 3.:  convex domains

class Convex_domains : public All_Domains {
 public:
  // Konstruktoren, n ist fuer Parameter n_uniform 
  Convex_domains(PointtypeD (*f)(D3vector),
                 double (*fd)(D3vector,dir_3D),
                 int n) {
    function = f;
    func_distance = fd;
    A = D3vector(0.0,0.0,0.0);
    H = 1.0;
    VecH.x = 1.0;
    VecH.y = 1.0;
    VecH.z = 1.0;
    n_uniform = n;
    is_periodic = false;
  }

  // Groesse eines uniformen Grundgitters, Angabe in Levels:
  int  Give_n_uniform() { return n_uniform; } ;
  // Kann berechnet werden ob Kante im Gebiet?
  bool calc_edge(D3vector V1, D3vector V2) { return true; }
  // ist Kante im Gebiet?
  bool edge(D3vector V1, D3vector V2);
  // Kann der Abstand zum Rand in Richting d berechnet werden?
  calc_dis calc_distance(D3vector V, dir_3D d) {return yes; }
  // Abstand zum Rand in Richting d
  double distance(D3vector V, dir_3D d);
  // ist Punkt im Gebiet? Diese Funktion sollte nur amAnfang verwendet werden.
  bool point_in_domain(D3vector V) {  return function(V); }

  // Periodoic object??
  bool Is_periodic() { return is_periodic; }

  // Beschreibung des Umgebenden Quadrates
  D3vector GiveA()    const { return A; }
  double   GiveH()    const { return H; }
  D3vector GiveVecH() const { return VecH; }
};

// convex domains with two parameters
class Convex_domains_2P : public All_Domains {
 public:
  // Konstruktoren, n ist fuer Parameter n_uniform 
  Convex_domains_2P(PointtypeD (*f)(D3vector,double,double),
                 double (*fd)(D3vector,dir_3D,double,double),
                 int n) {
    function_2P = f;
    func_distance_2P = fd;
    n_uniform = n;
    is_periodic = false;
  }

  // Groesse eines uniformen Grundgitters, Angabe in Levels:
  int  Give_n_uniform() { return n_uniform; } ;
  // Kann berechnet werden ob Kante im Gebiet?
  bool calc_edge(D3vector V1, D3vector V2) { return true; }
  // ist Kante im Gebiet?
  bool edge(D3vector V1, D3vector V2);
  // Kann der Abstand zum Rand in Richting d berechnet werden?
  calc_dis calc_distance(D3vector V, dir_3D d) {return yes; }
  // Abstand zum Rand in Richting d
  double distance(D3vector V, dir_3D d);
  // ist Punkt im Gebiet? Diese Funktion sollte nur amAnfang verwendet werden.
  bool point_in_domain(D3vector V) { 
    return function_2P(V,R_down,R_top);
  }

  // Periodoic object??
  bool Is_periodic() { return is_periodic; }

  // Beschreibung des Umgebenden Quadrates
  D3vector GiveA()    const { return A; }
  double   GiveH()    const { return H; }
  D3vector GiveVecH() const { return VecH; }
};

// convex domains with four parameters
class Convex_domains_4P : public All_Domains {
 public:
  // Konstruktoren, n ist fuer Parameter n_uniform 
  Convex_domains_4P(PointtypeD (*f)(D3vector,double,double,double,double),
                 double (*fd)(D3vector,dir_3D,double,double,double,double),
                 int n) {
    function_4P = f;
    func_distance_4P = fd;
    n_uniform = n;
    is_periodic = false;
  }

  // Groesse eines uniformen Grundgitters, Angabe in Levels:
  int  Give_n_uniform() { return n_uniform; } ;
  // Kann berechnet werden ob Kante im Gebiet?
  bool calc_edge(D3vector V1, D3vector V2) { return true; }
  // ist Kante im Gebiet?
  bool edge(D3vector V1, D3vector V2);
  // Kann der Abstand zum Rand in Richting d berechnet werden?
  calc_dis calc_distance(D3vector V, dir_3D d) {return yes; }
  // Abstand zum Rand in Richting d
  double distance(D3vector V, dir_3D d);
  // ist Punkt im Gebiet? Diese Funktion sollte nur amAnfang verwendet werden.
  bool point_in_domain(D3vector V) { 
    return function_4P(V,R_down,R_top,R_left,R_right);
  }

  // Periodoic object??
  bool Is_periodic() { return is_periodic; }

  // Beschreibung des Umgebenden Quadrates
  D3vector GiveA()    const { return A; }
  double   GiveH()    const { return H; }
  D3vector GiveVecH() const { return VecH; }
};

// 4.1:  union of domains

class Union_domains : public All_Domains {
 public:
  // constructors
  Union_domains(All_Domains* dom_a, All_Domains *dom_b) {
    D3vector A_a;
    D3vector A_b;
    D3vector VH_a;
    D3vector VH_b;
    domain_A = dom_a;
    domain_B = dom_b;
    A_a  = domain_A->GiveA();
    A_b  = domain_B->GiveA();
    VH_a = domain_A->GiveVecH();
    VH_b = domain_B->GiveVecH();
    A = D3vector(MIN(A_a.x,A_b.x),MIN(A_a.y,A_b.y),MIN(A_a.z,A_b.z));
    H = MAX(domain_A->GiveH(),domain_B->GiveH());
    VecH = D3vector(MAX(VH_a.x,VH_b.x),MAX(VH_a.y,VH_b.y),MAX(VH_a.z,VH_b.z));
    // x
    if(A_a.x < A_b.x) {
      H = MAX(A_b.x - A_a.x + domain_B->GiveH(),H);
      VecH.x = MAX(A_b.x - A_a.x + domain_B->GiveVecH().x,VecH.x);
    }
    if(A_a.x > A_b.x) {
      H = MAX(A_a.x - A_b.x + domain_A->GiveH(),H);
      VecH.x = MAX(A_a.x - A_b.x + domain_A->GiveVecH().x,VecH.x);
    }
    // y
    if(A_a.y < A_b.y) {
      H = MAX(A_b.y - A_a.y + domain_B->GiveH(),H);
      VecH.y = MAX(A_b.y - A_a.y + domain_B->GiveVecH().y,VecH.y);
    }
    if(A_a.y > A_b.y) {
      H = MAX(A_a.y - A_b.y + domain_A->GiveH(),H);
      VecH.y = MAX(A_a.y - A_b.y + domain_A->GiveVecH().y,VecH.y);
    }
    // z
    if(A_a.z < A_b.z) {
      H = MAX(A_b.z - A_a.z + domain_B->GiveH(),H);
      VecH.z = MAX(A_b.z - A_a.z + domain_B->GiveVecH().z,VecH.z);
    }
    if(A_a.z > A_b.z) {
      H = MAX(A_a.z - A_b.z + domain_A->GiveH(),H);
      VecH.z = MAX(A_a.z - A_b.z + domain_A->GiveVecH().z,VecH.z);
    }

    n_uniform = MAX(domain_A->Give_n_uniform(),domain_B->Give_n_uniform());
  }

  // Periodoic object??
  bool Is_periodic() { 
    return (domain_A->Is_periodic() || domain_B->Is_periodic());
  }
  
  // size of minimal uniform grid (in levels)
  int  Give_n_uniform() { return n_uniform; } ;
  // is it possible to calculate whether edge in domain?
  bool calc_edge(D3vector V1, D3vector V2) { return true; }
  // is edge in domain?
  bool edge(D3vector V1, D3vector V2);
  //  is it possible to calculate distance to the boundary?
  calc_dis calc_distance(D3vector V, dir_3D d) {return yes; }
  // distance to the boundary in direction d
  double distance(D3vector V, dir_3D d);
  // is point in domain?
  bool point_in_domain(D3vector V) {  
    return domain_A->point_in_domain(V)==interiorD ||
      domain_B->point_in_domain(V)==interiorD; }
  
  // describtion of the square, in which the domain is contained
  D3vector GiveA() const { return A; }
  double   GiveH() const { return H; }
  D3vector GiveVecH() const { return VecH; }
};



// 4.2:  cut of domains

class Cut_domains : public All_Domains {
 public:
  // constructors
  Cut_domains(All_Domains* dom_a, All_Domains *dom_b) {
    D3vector A_a;
    D3vector A_b;
    D3vector B_a;
    D3vector B_b;
    D3vector B;
    domain_A = dom_a;
    domain_B = dom_b;
    A_a = domain_A->GiveA();
    A_b = domain_B->GiveA();
    B_a = domain_A->GiveH() + A_a;
    B_b = domain_B->GiveH() + A_b;
    A = D3vector(MAX(A_a.x,A_b.x),MAX(A_a.y,A_b.y),MAX(A_a.z,A_b.z));
    B = D3vector(MIN(B_a.x,B_b.x),MIN(B_a.y,B_b.y),MIN(B_a.z,B_b.z));
    H = MAX(B.x-A.x,B.y-A.y,B.z-A.z);
    VecH = D3vector(B.x-A.x,B.y-A.y,B.z-A.z);
    //    H = MIN(domain_A->GiveH(),domain_B->GiveH());
    n_uniform = MAX(domain_A->Give_n_uniform(),domain_B->Give_n_uniform());
  }
  
  // size of minimal uniform grid (in levels)
  int  Give_n_uniform() { return n_uniform; } ;
  // is it possible to calculate whether edge in domain?
  bool calc_edge(D3vector V1, D3vector V2) { return true; }
  // is edge in domain?
  bool edge(D3vector V1, D3vector V2);
  //  is it possible to calculate distance to the boundary?
  calc_dis calc_distance(D3vector V, dir_3D d) {return yes; }
  // distance to the boundary in direction d
  double distance(D3vector V, dir_3D d);
  // is point in domain?
  bool point_in_domain(D3vector V) {  
    return domain_A->point_in_domain(V)==interiorD &&
      domain_B->point_in_domain(V)==interiorD; }
  
  // Periodoic object??
  bool Is_periodic() { 
    return (domain_A->Is_periodic() || domain_B->Is_periodic());
  }


  // describtion of the square, in which the domain is contained
  D3vector GiveA() const { return A; }
  double   GiveH() const { return H; }
  D3vector GiveVecH() const { return VecH; }
};


// 5.:  move domain

class Move_domain : public All_Domains {
 public:
  // constructor
  Move_domain(All_Domains *dom_a, D3vector Vmove)  {
    domain_A = dom_a;
    V_move = Vmove;
    A = domain_A->GiveA() + V_move;
    H = domain_A->GiveH();
    VecH = domain_A->GiveVecH();
    n_uniform = domain_A->Give_n_uniform();
  }

  // size of minimal uniform grid (in levels)
  int  Give_n_uniform() { return n_uniform; }
  // is it possible to calculate whether edge in domain?
  bool calc_edge(D3vector V1, D3vector V2) { return true; }
  // is edge in domain?
  bool edge(D3vector V1, D3vector V2);
  //  is it possible to calculate distance to the boundary?
  calc_dis calc_distance(D3vector V, dir_3D d) {return yes; }
  // distance to the boundary in direction d
  double distance(D3vector V, dir_3D d);
  // is point in domain? 
  bool point_in_domain(D3vector V) {  
    return domain_A->point_in_domain(V-V_move)==interiorD; }

  // Periodoic object??
  bool Is_periodic() { 
    return domain_A->Is_periodic();
  }


  // describtion of the square, in which the domain is contained
  D3vector GiveA() const { return A; }
  double   GiveH() const { return H; }
  D3vector GiveVecH() const { return VecH; }
};

// 6.:  scale domain

class Scale_domain : public All_Domains {
 public:
  // constructor
  Scale_domain(All_Domains *dom_a, D3vector Vscale)  {
    domain_A = dom_a;
    V_scale.x = ABS(Vscale.x);
    V_scale.y = ABS(Vscale.y);
    V_scale.z = ABS(Vscale.z);
    A = domain_A->GiveA() * Vscale;
    H = domain_A->GiveH() * MAX(V_scale.x,V_scale.y,V_scale.z);
    VecH = domain_A->GiveVecH() * D3vector(V_scale.x,V_scale.y,V_scale.z);
    n_uniform = domain_A->Give_n_uniform();
  }

  // size of minimal uniform grid (in levels)
  int  Give_n_uniform() { return n_uniform; }
  // is it possible to calculate whether edge in domain?
  bool calc_edge(D3vector V1, D3vector V2) { return true; }
  // is edge in domain?
  bool edge(D3vector V1, D3vector V2);
  //  is it possible to calculate distance to the boundary?
  calc_dis calc_distance(D3vector V, dir_3D d) {return yes; }
  // distance to the boundary in direction d
  double distance(D3vector V, dir_3D d);
 // is point in domain? 
  bool point_in_domain(D3vector V) {  
    return domain_A->point_in_domain(V/V_scale)==interiorD; }

  // Periodoic object??
  bool Is_periodic() { 
    return domain_A->Is_periodic();
  }

  // describtion of the square, in which the domain is contained
  D3vector GiveA() const { return A; }
  double   GiveH() const { return H; }
  D3vector GiveVecH() const { return VecH; }
};

// 7.:  operators
// operator +
inline Domain
operator+(All_Domains& a, D3vector V)
{
  return Domain(new Move_domain(&a,V));
}

inline Domain
operator+(D3vector V, All_Domains& a)
{
  return Domain(new Move_domain(&a,V));
}

inline Domain
operator+(const Domain& a, D3vector V)
{
  return Domain(new Move_domain(new Domain(a),V));
}

inline Domain
operator+(D3vector V, const Domain& a)
{
  return Domain(new Move_domain(new Domain(a),V));
}
// operator *
inline Domain
operator*(All_Domains& a, const D3vector V)
{
  return Domain(new Scale_domain(&a,V));
}

inline Domain
operator*(const Domain& a, const D3vector V)
{
  return Domain(new Scale_domain(new Domain(a),V));
}

inline Domain
operator*(All_Domains& a, double s)
{
  return Domain(new Scale_domain(&a,D3vector(s,s,s)));
}

inline Domain
operator*(double s, All_Domains& a)
{
  return Domain(new Scale_domain(&a,D3vector(s,s,s)));
}

inline Domain
operator*(const Domain& a, double s)
{
  return Domain(new Scale_domain(new Domain(a),D3vector(s,s,s)));
}

inline Domain
operator*(double s, const Domain& a)
{
  return Domain(new Scale_domain(new Domain(a),D3vector(s,s,s)));
}

// operator ||
inline Domain
operator||(All_Domains& a, const Domain& b)
{
  return Domain(new Union_domains(&a,new Domain(b)));
}

inline Domain
operator||(All_Domains& a, All_Domains& b)
{
  return Domain(new Union_domains(&a,&b));
}

inline Domain
operator||(const Domain& a, All_Domains& b)
{
  return Domain(new Union_domains(new Domain(a),&b));
}

inline Domain
operator||(const Domain& a, const Domain& b)
{
  return Domain(new Union_domains(new Domain(a),new Domain(b)));
}

// operator &&
inline Domain
operator&&(All_Domains& a, const Domain& b)
{
  return Domain(new Cut_domains(&a,new Domain(b)));
}

inline Domain
operator&&(All_Domains& a, All_Domains& b)
{
  return Domain(new Cut_domains(&a,&b));
}

inline Domain
operator&&(const Domain& a, All_Domains& b)
{
  return Domain(new Cut_domains(new Domain(a),&b));
}

inline Domain
operator&&(const Domain& a, const Domain& b)
{
  return Domain(new Cut_domains(new Domain(a),new Domain(b)));
}


// some elementary member functions
// --------------------------------------
// Convex_domains
inline double Convex_domains::distance(D3vector V, dir_3D d) {
  return func_distance(V,d);
}

inline bool Convex_domains::edge(D3vector V1, D3vector V2) {
  if(function(V1) == exteriorD || function(V2) == exteriorD)
    return false;
  else return true;
}

// Convex_domains
inline double Convex_domains_2P::distance(D3vector V, dir_3D d) {
  return func_distance_2P(V,d,R_down,R_top);
}
inline double Convex_domains_4P::distance(D3vector V, dir_3D d) {
  return func_distance_4P(V,d,R_down,R_top,R_left,R_right);
}

inline bool Convex_domains_2P::edge(D3vector V1, D3vector V2) {
  if(function_2P(V1,R_down,R_top) == exteriorD ||
     function_2P(V2,R_down,R_top) == exteriorD)
    return false;
  else return true;
}

inline bool Convex_domains_4P::edge(D3vector V1, D3vector V2) {
  if(function_4P(V1,R_down,R_top,R_left,R_right) == exteriorD ||
     function_4P(V2,R_down,R_top,R_left,R_right) == exteriorD)
    return false;
  else return true;
}

// Move_domains
inline double Move_domain::distance(D3vector V, dir_3D d) {
  return domain_A->distance(V-V_move,d);
}

inline bool Move_domain::edge(D3vector V1, D3vector V2) {
  return domain_A->edge(V1-V_move,V2-V_move);
}

// Scale_domains
inline double Scale_domain::distance(D3vector V, dir_3D d) {
  if(d<Sdir) return domain_A->distance(V/V_scale,d) * V_scale.x;
  if(d>Ndir) return domain_A->distance(V/V_scale,d) * V_scale.z;
  return domain_A->distance(V/V_scale,d) * V_scale.y;
}

inline bool Scale_domain::edge(D3vector V1, D3vector V2) {
  return domain_A->edge(V1/V_scale,V2/V_scale);
}

// Union_domains
inline bool Union_domains::edge(D3vector V1, D3vector V2) {
  double eps, dis;
  if(domain_A->edge(V1-V_move,V2-V_move)) return true;
  if(domain_B->edge(V1-V_move,V2-V_move)) return true;
  eps = H * domain_eps;
  // E - W
  dis = V1.x - V2.x;
  if(dis > eps) {  //  V2.x ---- V1.x
    if(dis < domain_A->distance(V2,Edir) + domain_B->distance(V1,Wdir))
      return true;
    if(dis < domain_B->distance(V2,Edir) + domain_A->distance(V1,Wdir))
      return true;
    return false;
  }
  dis = V2.x - V1.x;
  if(dis > eps) {  //  V1.x ---- V2.x
    if(dis < domain_A->distance(V1,Edir) + domain_B->distance(V2,Wdir))
      return true;
    if(dis < domain_B->distance(V1,Edir) + domain_A->distance(V2,Wdir))
      return true;
    return false;
  }
  // S - N
  dis = V1.y - V2.y;
  if(dis > eps) {  //  V2.y ---- V1.y
    if(dis < domain_A->distance(V2,Ndir) + domain_B->distance(V1,Sdir))
      return true;
    if(dis < domain_B->distance(V2,Ndir) + domain_A->distance(V1,Sdir))
      return true;
    return false;
  }
  dis = V2.y - V1.y;
  if(dis > eps) {  //  V1.y ---- V2.y
    if(dis < domain_A->distance(V1,Ndir) + domain_B->distance(V2,Sdir))
      return true;
    if(dis < domain_B->distance(V1,Ndir) + domain_A->distance(V2,Sdir))
      return true;
    return false;
  }
  // D - T
  dis = V1.z - V2.z;
  if(dis > eps) {  //  V2.z ---- V1.z
    if(dis < domain_A->distance(V2,Tdir) + domain_B->distance(V1,Ddir))
      return true;
    if(dis < domain_B->distance(V2,Tdir) + domain_A->distance(V1,Ddir))
      return true;
    return false;
  }
  dis = V2.z - V1.z;
  //  if(dis > eps) {  //  V1.z ---- V2.z
  if(dis < domain_A->distance(V1,Tdir) + domain_B->distance(V2,Ddir))
    return true;
  if(dis < domain_B->distance(V1,Tdir) + domain_A->distance(V2,Ddir))
    return true;
  return false;
  // }
}

inline double Union_domains::distance(D3vector V, dir_3D d) {
  double dA, dB, eps;
  eps = H * domain_eps;
  dA = domain_A->distance(V,d);
  dB = domain_B->distance(V,d);
  if(!(dA>0.0) && dB>0.0)
    return MAX(domain_A->distance(V+(dB-eps)*Unit_vector(d),d)+dB+eps,dB);
  if(dA>0.0 && !(dB>0.0))
    return MAX(dA,domain_B->distance(V+(dA-eps)*Unit_vector(d),d)+dA+eps);
  return MAX(domain_A->distance(V,d),
             domain_B->distance(V,d));
}



// Cut_domains
inline bool Cut_domains::edge(D3vector V1, D3vector V2) {
  if(domain_A->edge(V1-V_move,V2-V_move) &&
     domain_B->edge(V1-V_move,V2-V_move)) return true;
  else return false;
}

inline double Cut_domains::distance(D3vector V, dir_3D d) {
  double dA, dB;
  dA = domain_A->distance(V,d);
  dB = domain_B->distance(V,d);
  return MIN(dA,dB);
}



// Remark:
// the functions
//        int Give_n_uniform() =0;
// are not implemented in an optimal way.
// sometimes they give back too small values
// this is not a problem, since this function is
// used only for error messages!!

// 9.:  special domain fuer Andreas  


class Special_domain : public Convex_domains {
 public:
  Special_domain(double (*dis_Specdom)(D3vector V, dir_3D d), PointtypeD (*Poi_Specdom)(D3vector V), 
                 double transverseScale);
};

inline Special_domain ::Special_domain (double (*dis_Specdom)(D3vector V, dir_3D d), PointtypeD (*Poi_Specdom)(D3vector V),
                                        double transverseScale) : 
  Convex_domains(Poi_Specdom,dis_Specdom,1) {
  A = D3vector(-transverseScale,-transverseScale,0.0);
  H = 2.0*transverseScale;

  VecH = D3vector(2.0*transverseScale,2.0*transverseScale,1.0);
}


#endif
    

---