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

// ------------------------------------------------------------
// math_lib.h
//
// ------------------------------------------------------------

#ifndef MATHLIB_H_
#define MATHLIB_H_
       
#ifndef M_PI
#define M_PI 3.1415927
#endif

#ifdef COMP_GNUOLD

#else
 #include <iostream>
 #include <fstream>           
 using std::ofstream;
 using std::cout;
 using std::endl;           
#endif 
         

// 1. next prim
// 2. MAX, MIN
// 3. 3D vector class
// 4. geometric operators for 3D vectors
// 5. others



// 1. next prim 


// Find next prim
int Find_next_prime(int zahl);


// trigonometric functions

inline double arccos(const double y) {
   return atan2(sqrt(1-y*y),y);
   // if(y>=0) return atan(sqrt(1-y*y)/y);
   // else    return atan(sqrt(1-y*y)/y) + M_PI;
}

inline double my_tan(const double y) {
  return tan(y);
}

inline double my_atan(const double y) {
  return atan(y);
}

inline double my_sqrt(const double y) {
  return sqrt(y);
}

// MAX, MIN


/* Maximum */
inline double MAX(double a, double b) {
  if(a < b) return b;
  else return a;
}

inline double MAX(double a, double b, double c) {
  return MAX(MAX(a,b),c);
}

inline int MAX(int a, int b) {
  if(a < b) return b;
  else return a;
}

inline int MAX(int a, int b, int c) {
  static int m;
  if(a < b) m = b;
  else m = a;
  if(c < m) return m;
  else return c;
}

/* Minimum */
inline double MIN(double a, double b) {
  if(a < b) return a;
  else return b;
}

/* Absolute value */
inline double ABS(double a) {
  if(a < 0.0) return -a;
  else return a;
}

/* Absolute value */
inline int ABS(int a) {
  if(a < 0) return -a;
  else return a;
}


// 3. a simple 3D vector class

class D3vector {
 public:
  double x,y,z;
  D3vector(double cx, double cy, double cz) : x(cx), y(cy), z(cz) {}; 
  D3vector(double c) : x(c), y(c), z(c) {}; 
  D3vector() : x(0), y(0), z(0) {}; 
  void Print();
  void Print(ofstream *Datei);
  void operator=(const D3vector& v) { x=v.x; y=v.y; z=v.z; }
  double operator[](int i) const {
    if(i==0) return x;
    if(i==1) return y;
    return z;
  }
};   

inline D3vector MAX(D3vector V1, D3vector V2) {
  return D3vector(MAX(V1.x,V2.x),MAX(V1.y,V2.y),MAX(V1.z,V2.z));
}

inline D3vector MIN(D3vector V1, D3vector V2) {
  return D3vector(MIN(V1.x,V2.x),MIN(V1.y,V2.y),MIN(V1.z,V2.z));
}

inline D3vector operator+(const D3vector& v,const D3vector& w) {
  return D3vector(v.x+w.x,v.y+w.y,v.z+w.z);
};

inline D3vector operator*(const D3vector& v,const D3vector& w) {
  return D3vector(v.x*w.x,v.y*w.y,v.z*w.z);
};

inline D3vector operator/(const D3vector& v,const D3vector& w) {
  return D3vector(v.x/w.x,v.y/w.y,v.z/w.z);
};

inline D3vector operator-(const D3vector& v,const D3vector& w) {
  return D3vector(v.x-w.x,v.y-w.y,v.z-w.z);
};

inline D3vector operator/(const D3vector& v,const double f) {
  return D3vector(v.x/f,v.y/f,v.z/f);
};

inline D3vector operator*(const D3vector& v,const double f) {
  return D3vector(v.x*f,v.y*f,v.z*f);
};

inline D3vector operator*(const double f, const D3vector& v) {
  return D3vector(v.x*f,v.y*f,v.z*f);
};

inline bool operator<(const D3vector& v,const D3vector& w) {
  return (v.x<w.x && v.y<w.y && v.z<w.z);
};
inline bool operator==(const D3vector& v,const D3vector& w) {
  return (v.x==w.x && v.y==w.y && v.z==w.z);
};


// 4. geometric operators for 3D vectors

inline double L_infty(const D3vector& v) {
  return MAX(ABS(v.x),ABS(v.y),ABS(v.z));
}

inline double norm(const D3vector& v) {
  return my_sqrt(v.x*v.x+v.y*v.y+v.z*v.z);
}

// scalar product
inline double product(const D3vector& v, const D3vector& w) {
  return v.x*w.x+v.y*w.y+v.z*w.z;
}

inline D3vector cross_product(const D3vector& v, const D3vector& w) {
  return D3vector(v.y*w.z-v.z*w.y,
                  v.z*w.x-v.x*w.z,
                  v.x*w.y-v.y*w.x);
};

inline D3vector normal_vector_of_triangle(const D3vector& va,
                                          const D3vector& vb,
                                          const D3vector& vc) {
  D3vector z;
  z = cross_product(va-vb,vc-vb);
  return z / norm(z);
}


inline double angle_between_vectors(const D3vector& va,
                                    const D3vector& vb) {
  return arccos(product(va,vb) / (norm(va)*norm(vb))) / M_PI * 180 ;
}



inline double max_interior_angel_of_triangle(const D3vector& va,
                                             const D3vector& vb,
                                             const D3vector& vc) {
  return MAX(angle_between_vectors(vc-va,vb-va),
             angle_between_vectors(va-vb,vc-vb),
             angle_between_vectors(vb-vc,va-vc));
}

// va, vb, vc is one face
// vb, vd, vc is one face
inline double angle_between_faces(const D3vector& va,
                                  const D3vector& vb,
                                  const D3vector& vc,
                                  const D3vector& vd) {
  return 180 - angle_between_vectors(normal_vector_of_triangle(va,vb,vc),
                                     normal_vector_of_triangle(vb,vd,vc));
}

double calc_maximal_edge_angle(const D3vector& va,
                               const D3vector& vb,
                               const D3vector& vc,
                               const D3vector& vd);

double calc_maximal_face_angle(const D3vector& va,
                               const D3vector& vb,
                               const D3vector& vc,
                               const D3vector& vd);

// volumne of tetrahedra

#define x0 V0.x
#define y0 V0.y
#define z0 V0.z

#define x1 V1.x
#define y1 V1.y
#define z1 V1.z

#define x2 V2.x
#define y2 V2.y
#define z2 V2.z

#define x3 V3.x
#define y3 V3.y
#define z3 V3.z

inline double Calc_det(const D3vector& V0, const D3vector& V1,
                       const D3vector& V2, const D3vector& V3) {
  return (x1*y2*z0 - 
          x1*y3*z0 - x0*y2*z1 + x0*y3*z1 - x1*y0*z2 + x0*y1*z2 - x0*y3*z2 + 
          x1*y3*z2 + x3*(y1*z0 - y2*z0 - y0*z1 + y2*z1 + y0*z2 - y1*z2) + 
          x1*y0*z3 - x0*y1*z3 + x0*y2*z3 - x1*y2*z3 + 
          x2*(y3*z0 + y0*z1 - y3*z1 - y0*z3 + y1*(-z0 + z3)));
}

#undef x0
#undef y0
#undef z0

#undef x1
#undef y1
#undef z1

#undef x2
#undef y2
#undef z2

#undef x3
#undef y3
#undef z3

inline double Calc_bary_weights(const D3vector& V0, const D3vector& V1,
                                const D3vector& V2, const D3vector& V3,
                                const D3vector& V, int i) {
  if(i==0) return Calc_det(V,V1,V2,V3);
  if(i==1) return Calc_det(V0,V,V2,V3);
  if(i==2) return Calc_det(V0,V1,V,V3);
  if(i==3) return Calc_det(V0,V1,V2,V);
  cout << " error in Calc_bary_weights! " << endl;
  return (double) -1.0;
}



// 5. others


/* Test ob Zahl gerade */
inline int gerade(int N) {
  return (~N)&1;
}

// Berechnung der Zweierpotenz 
inline int Zweipotenz (int k) {
  return 1<<k;
}


// Implementierung einiger Memberfunktionen


// D3vector 
// ----------

inline void D3vector::Print() {
  cout << "Coordinate: " << x << ", " << y << ", " << z << ";";
};

inline void D3vector::Print(ofstream *Datei) {
  *Datei  << x << " " << y << " " << z;
};


     

#endif
   
  

---