Helpers: static method -> namespace methods

extensions/GeometryTools.cc

@@ -8,39 +8,6 @@ namespace meshconv {
 //                                      ##  General  ##
 //                                      ###############
 
-//----------------------------------------< coordinate_extrema_triangle >
-//returns in 'small_coord' the smallest coordinate of 'sm' in the given dimension and in 'big_coord' the
-//biggest coordinate in the given dimension. Note: 'sm' has to be a triangle in a 2d or 3d world.
-// void coordinate_extrema_triangle(mesh &m, element &sm, long dimension, double &small_coord,
-//                                                                           double &big_coord){
-//   if(m.vertices[sm.v[0]][dimension] <= m.vertices[sm.v[1]][dimension]){
-//     if(m.vertices[sm.v[0]][dimension] <= m.vertices[sm.v[2]][dimension]){
-//       small_coord = m.vertices[sm.v[0]][dimension];
-//       if(m.vertices[sm.v[1]][dimension] >= m.vertices[sm.v[2]][dimension])
-//         big_coord = m.vertices[sm.v[1]][dimension];
-//       else
-//         big_coord = m.vertices[sm.v[2]][dimension];
-//     }
-//     else{
-//       small_coord = m.vertices[sm.v[2]][dimension];
-//       big_coord = m.vertices[sm.v[1]][dimension];
-//     }
-//   }
-//   else{
-//     if(m.vertices[sm.v[1]][dimension] <= m.vertices[sm.v[2]][dimension]){
-//       small_coord = m.vertices[sm.v[1]][dimension];
-//       if(m.vertices[sm.v[0]][dimension] >= m.vertices[sm.v[2]][dimension])
-//         big_coord = m.vertices[sm.v[0]][dimension];
-//       else
-//         big_coord = m.vertices[sm.v[2]][dimension];
-//     }
-//     else{
-//       small_coord = m.vertices[sm.v[2]][dimension];
-//       big_coord = m.vertices[sm.v[0]][dimension];
-//     }
-//   }
-// }
-
 //----------------------------------------< solve_determinant4 >
 double solve_determinant4(double x11, double x12, double x13, double x14,
                           double x21, double x22, double x23, double x24,
@@ -447,6 +414,11 @@ bool point_in_polygon(double point[], const std::vector<AMDiS::WorldVector<doubl
   return inside;
 }
 
+//----------------------------------------< triangle_area_2d >
+double triangle_area_2d(double tri0[], double tri1[], double tri2[]){
+  return 0.5 * abs((tri1[0] - tri0[0]) * (tri2[1] - tri0[1]) - (tri2[0] - tri0[0]) * (tri1[1] - tri0[1]));
+}
+
 //                                      #################
 //                                      ##  3D Worlds  ##
 //                                      #################
@@ -1500,4 +1472,14 @@ double distance_point_triangle_3d(double point[], double tri0[], double tri1[],
   return distance;
 }
 
+
+double volume_tetraeder(double tet0[], double tet1[], double tet2[], double tet3[]){
+  mtl::dense2d<double> A(3,3);
+  A(0,0) = tet1[0]-tet0[0]; A(0,1) = tet2[0]-tet0[0]; A(0,2) = tet3[0]-tet0[0];
+  A(1,0) = tet1[1]-tet0[1]; A(1,1) = tet2[1]-tet0[1]; A(1,2) = tet3[1]-tet0[1];
+  A(2,0) = tet1[2]-tet0[2]; A(2,1) = tet2[2]-tet0[2]; A(2,2) = tet3[2]-tet0[2];
+    
+  return abs(Helpers::determinant(A))/6.0;
+}
+
 }

extensions/GeometryTools.h

@@ -93,6 +93,9 @@ double distance_point_triangle_2d(double point[], double tri0[], double tri1[],
 
 bool point_in_polygon(double point[], const std::vector<AMDiS::WorldVector<double> > &vertices);
 
+//----------------------------------------< triangle_area_2d >
+double triangle_area_2d(double tri0[], double tri1[], double tri2[]);
+
 //                                      #################
 //                                      ##  3D Worlds  ##
 //                                      #################
@@ -206,45 +209,13 @@ double distance_point_box_3d(double point[], double min_corner[], double max_cor
 //calculates the distance between a given 'point' and the triangle given by 'tri0', 'tri1', 'tri2'.
 double distance_point_triangle_3d(double point[], double tri0[], double tri1[], double tri2[]);
 
-
 template<int dow>
-void coordinate_transform(double x[], double shift[], double scale[], double alpha=0.0, double beta=0.0)
-{
-  mtl::dense2D<double> Rx(dow,dow),Ry(dow,dow),T(dow,dow),S(dow,dow);
-  Rx = 1.0;
-  Ry = 1.0;
-  S = 1.0;
-  if (dow == 3) {
-    Rx(1,1) = cos(alpha);  Rx(1,2) = -sin(alpha);
-    Rx(2,1) = sin(alpha);  Rx(2,2) = cos(alpha);
-
-    Ry(0,0) = cos(beta);  Ry(0,2) = sin(beta);
-    Ry(2,0) = -sin(beta); Ry(2,2) = cos(beta);
-
-    S(0,0) = scale[0];  S(1,1) = scale[1];  S(2,2) = scale[2];
-  } else if (dow == 2) {
-    Rx(0,0) = cos(alpha);  Rx(0,1) = -sin(alpha);
-    Rx(1,0) = sin(alpha);  Rx(1,1) = cos(alpha);
-
-    Ry = 1.0;
-
-    S(0,0) = scale[0];  S(1,1) = scale[1];
-  } else {
-    S(0,0) = scale[0];
-  }
-
-  T = Ry*Rx*S;
-
-  mtl::dense_vector<double> coords(dow);
-  for (int i = 0; i < dow; i++)
-    coords[i] = x[i];
+void coordinate_transform(double x[], double shift[], double scale[], double alpha = 0.0, double beta = 0.0);
 
-  coords = T*coords;
+double volume_tetraeder(double tet0[], double tet1[], double tet2[], double tet3[]);
 
-  for (int i = 0; i < dow; i++)
-    x[i] = coords[i]+shift[i];
-}
+} // end namespace meshconv
 
-}
+#include "GeometryTools.hh"
 
 #endif // GEOMETRY_TOOLBOX_H

extensions/Helpers.h

@@ -22,44 +22,48 @@
 #include <stdarg.h>
 
 #include "VectorOperations.h"
+#include "Views.h"
 
 #define TEST_WARNING(test) if ((test));else WARNING
 
 using namespace std;
 using namespace AMDiS; 
 
-static long random_seed_initial_value=0;
+static long random_seed_initial_value = 0;
 
-class Helpers {
-public:
+namespace Helpers {
 
   /// math routines
   /// ==============
   
-  static double cint(double x){
+  inline double cint(double x) {
     double intpart; 
-    if (modf(x,&intpart)>=.5)
-      return x>=0?ceil(x):floor(x);
+    if (modf(x, &intpart) >= 0.5)
+      return x >=0 ? ceil(x) : floor(x);
     else
-      return x<0?ceil(x):floor(x);
+      return x < 0 ? ceil(x) : floor(x);
   }
   
-  static double round(double r, double places){
-    double off= pow(10.0, places);
+  inline double round (double r, double places) {
+    double off = pow(10.0, places);
     return cint(r*off)/off;
   }
   
-  static bool isNumber(double val) { return !isnan(val) && !isinf(val); };
-  static double signum(double val, double tol=1.e-10) { return (val>tol?1.0:(val<-tol?-1.0:0.0)); }
-  static int signum(int val) { return (val>0?1:(val<0?-1:0)); }
-  static int signum(bool val) { return (val?1:-1); }
-  static int toInt(bool val) { return (val?1:0); }
+  inline double signum(double val, double tol = DBL_TOL) { return val > tol ? 1.0 : (val < -tol ? -1.0 : 0.0); }
+  inline int signum(int val) { return val > 0 ? 1 : (val < 0 ? -1 : 0); }
+  inline int signum(bool val) { return val ? 1 : -1; }
+  inline int toInt(bool val) { return val ? 1 : 0; }
+  inline bool isNumber(double val) { return !isnan(val) && !isinf(val); }
+  inline bool toBool(double val, double tol = DBL_TOL) { return val > tol || val < -tol ? true : false; }
+  inline bool isPositiv(double val, double tol = DBL_TOL) { return val > -tol; }
+  inline bool isStrictPositiv(double val, double tol = DBL_TOL) { return val > tol; }
 
   /// routines for conversion to string
   /// =================================
   
   template<typename T>
-  static std::string toString(const T &value, ios_base::fmtflags flag= ios_base::scientific)
+  inline std::string toString(const T &value, 
+			      ios_base::fmtflags flag = ios_base::scientific)
   {
     std::stringstream ss;
     ss.setf(flag);
@@ -67,15 +71,17 @@ public:
     if (!(ss<<value))
       throw(std::runtime_error("toString() not possible"));
     return ss.str();
-  };
+  }
   
-  static std::string toString(const int &value)
+  inline std::string toString(const int &value)
   {
     return boost::lexical_cast<std::string>(value);
-  };
+  }
   
   template<typename T>
-  static std::string toString(const WorldVector<T> &value, bool brackets= true, ios_base::fmtflags flag= ios_base::scientific)
+  inline std::string toString(const WorldVector<T> &value, 
+			      bool brackets = true, 
+			      ios_base::fmtflags flag = ios_base::scientific)
   {
     std::stringstream ss;
     if (brackets) ss<<"[";
@@ -86,10 +92,12 @@ public:
     
     if (brackets) ss<<"]";
     return ss.str();
-  };
+  }
   
   template<typename T>
-  static std::string toString(const std::vector<T> &vec, bool brackets= true, ios_base::fmtflags flag= ios_base::scientific)
+  inline std::string toString(const std::vector<T> &vec, 
+			      bool brackets = true, 
+			      ios_base::fmtflags flag = ios_base::scientific)
   {
     std::stringstream ss;
     if (brackets) ss<<"[";
@@ -101,18 +109,18 @@ public:
     }
     if (brackets) ss<<"]";
     return ss.str();
-  };
+  }
   
   template<class T>
-  static T fromString(const std::string& s)
+  inline T fromString(const std::string& s)
   {
     std::istringstream stream(s);
     T t;
     stream >> t;
     return t;
-  };
+  }
 
-  static std::string fillString(int length, char c, int numValues, ...)
+  inline std::string fillString(int length, char c, int numValues, ...)
   {
     va_list values;
     int value;
@@ -132,10 +140,10 @@ public:
     va_end(values);
 
     return std::string(len, c);
-  };
+  }
 
   // process printable string of memory
-  static std::string memoryToString(double mem, int startIdx=0) {
+  inline std::string memoryToString(double mem, int startIdx = 0) {
     int idx = startIdx;
     double mem_ = mem;
     while (mem_/1024.0 > 1.0) {
@@ -144,63 +152,53 @@ public:
     }
     std::string result = toString(mem_)+" "+(idx==0?"B":(idx==1?"KB":(idx==2?"MB":"GB")));
     return result;
-  };  
+  }
   
   
   /// some mesh routines
   /// ===========================
   
-  static void transformMesh(Mesh *mesh, WorldVector<double> scale, WorldVector<double> shift, WorldVector<double> rotate);
+  void transformMesh(Mesh *mesh, WorldVector<double> scale, WorldVector<double> shift, WorldVector<double> rotate);
 
-  static void scaleMesh(std::vector<Mesh*> meshes, WorldVector<double> scale);
+  void scaleMesh(std::vector<Mesh*> meshes, WorldVector<double> scale);
   // scale and shift by different values in all directions
-  static void scaleMesh(Mesh *mesh, WorldVector<double> shift, WorldVector<double> scale);
+  void scaleMesh(Mesh *mesh, WorldVector<double> shift, WorldVector<double> scale);
   // scale by different values in all directions
-  static void scaleMesh(Mesh *mesh, WorldVector<double> scale);
+  void scaleMesh(Mesh *mesh, WorldVector<double> scale);
   // scale and shift by the same values in all directions
-  static void scaleMesh(Mesh *mesh, double shift=0.0, double scale=1.0);
+  void scaleMesh(Mesh *mesh, double shift=0.0, double scale=1.0);
   
   /// calculate the dimension of a mesh
-  static WorldVector<double> getMeshDimension(Mesh *mesh);
-  static void getMeshDimension(Mesh *mesh, WorldVector<double> &min_corner, WorldVector<double> &max_corner);
+  WorldVector<double> getMeshDimension(Mesh *mesh);
+  void getMeshDimension(Mesh *mesh, WorldVector<double> &min_corner, WorldVector<double> &max_corner);
 
   /// read DOFVector from AMDiS .dat-files
-  static void read_dof_vector(const std::string file, DOFVector<double> *dofvec, long size);
-
-  /// copy DOFVectors that live on different meshes
-  static void copyDOF(std::vector<DOFVector<double>*> dof1, std::vector<DOFVector<double>*> dof2, WorldVector<double> shift, double fillValue = 0.0);
-  
-  static void copyDOF(DOFVector<double>* dof1, DOFVector<double>* dof2) {
-    WorldVector<double> shift; shift.set(0.0);
-    std::vector<DOFVector<double>*> dof1Vec; dof1Vec.push_back(dof1);
-    std::vector<DOFVector<double>*> dof2Vec; dof2Vec.push_back(dof2);
-    copyDOF(dof1Vec, dof2Vec, shift);
-  };
+  void read_dof_vector(const std::string file, DOFVector<double> *dofvec, long size);
 
   /// some linear algebra methods
   /// ===========================
 
-  static double determinant(mtl::dense2D<double> &A) // only for 3x3 - matrices
+  inline double determinant(mtl::dense2D<double> &A) // only for 3x3 - matrices
   {
     if(num_rows(A)==3 && num_cols(A)==3) {
       double det = A(0,0)*A(1,1)*A(2,2)+A(0,1)*A(1,2)*A(2,0)+A(0,2)*A(1,0)*A(2,1);
       return det-(A(0,2)*A(1,1)*A(2,0)+A(0,1)*A(1,0)*A(2,2)+A(0,0)*A(1,2)*A(2,1));
     }
     return 1.0;
-  };
+  }
   
-  static double determinant(WorldMatrix<double> &A) // only for 3x3 - matrices
+  inline double determinant(WorldMatrix<double> &A) // only for 3x3 - matrices
   {
     if(A.getNumRows()==3 && A.getNumCols()==3) {
       double det = A[0][0]*A[1][1]*A[2][2]+A[0][1]*A[1][2]*A[2][0]+A[0][2]*A[1][0]*A[2][1];
       return det-(A[0][2]*A[1][1]*A[2][0]+A[0][1]*A[1][0]*A[2][2]+A[0][0]*A[1][2]*A[2][1]);
     }
     return 1.0;
-  };
+  }
   
   /// inverse power method to find minimal eigenvalue of matrix A and appropriate eigenvector x, using m iterations
   template <class LinearSolver, class Matrix, class EigenVector>
-  static double inverse_iteration(LinearSolver& solver, const Matrix& A, EigenVector& x, int m)
+  inline double inverse_iteration(LinearSolver& solver, const Matrix& A, EigenVector& x, int m)
   { FUNCNAME("Helpers::inverse_iteration()");
 
     EigenVector y( size(x) ), res( size(x) );
@@ -212,20 +210,21 @@ public:
     solver.setMultipleRhs(true);
     for (int i = 0; i < m; ++i) {
       solver.solveSystem(A, y, x ); // solve Ay=x --> y
-      lambda_old=lambda;
-      lambda = dot(y,x)/dot(y,y);
-      relErr = abs((lambda-lambda_old)/lambda);
-      if(relErr<1.e-5) break;
+      lambda_old = lambda;
+      lambda = dot(y,x) / dot(y,y);
+      relErr = abs((lambda - lambda_old) / lambda);
+      if (relErr < 1.e-5)
+	break;
       x = y / two_norm(y);
     }
     solver.setMultipleRhs(false);
 
-    return lambda ;
-  };
+    return lambda;
+  }
 
   /// calculate approximation to condition number of matrix A using the OEMSolver solver
   template <class LinearSolver, class Matrix>
-  static double condition(LinearSolver& solver, const Matrix& A, int m=10)
+  inline double condition (LinearSolver& solver, const Matrix& A, int m=10)
   { FUNCNAME("Helpers::condition()");
 
     mtl::dense_vector<double>   x(num_rows(A)), y(num_rows(A));
@@ -236,22 +235,40 @@ public:
     x /= two_norm(x);
     for (int i = 0; i < m; ++i) {
       y = A * x;
-      lambda_old=lambda;
+      lambda_old = lambda;
       lambda = mtl::dot(y,x);
       relErr = abs((lambda-lambda_old)/lambda);
-      if(relErr<1.e-5) break;
+      if (relErr < 1.e-5)
+	break;
       x = y / two_norm(y) ;
     }
     result1 = std::abs(lambda / inverse_iteration(solver, A, x, m));
     return result1;
-  };
+  }
   
   /// interpolate values of DOFVector along line from (x1,y1) -> (x2,y1) with nPoints points
   template<typename Container>
-  static void interpolOverLine(const Container &vec, 
+  inline void interpolOverLine(const Container &vec, 
 			       double x1, double y1, double x2, double y2, 
 			       int nPoints,
-			       std::vector<WorldVector<double> > &x, std::vector<double> &y);
+			       std::vector<WorldVector<double> > &x, std::vector<double> &y)
+  { FUNCNAME("Helpers::interpolOverLine()");
+
+    WorldVector<double> p;
+
+    if (nPoints <= 1)
+      throw(std::runtime_error("Zu wenig Punkte fuer eine Interpolation!"));
+    
+    for (int i = 0; i < nPoints; ++i) {
+      double lambda = static_cast<double>(i)/static_cast<double>(nPoints - 1.0);
+      p[0] = lambda*x2 + (1.0-lambda)*x1;
+      p[1] = lambda*y2 + (1.0-lambda)*y1;
+      
+      double value = evalAtPoint(vec, p);
+      x.push_back(p);
+      y.push_back(value);
+    }
+  }
 
   /// calculate maxima of DOFVector along line, using interpolOverLine
   static void calcMaxOnXAxis(DOFVector<double> *rho, std::vector<std::pair<WorldVector<double>, double> > &maxima);
@@ -269,16 +286,16 @@ public:
   /// misc routines
   /// =============
   
-  static void plot(std::vector<double> &values, int numRows=7, int numCols=20, std::string symbol="*");
+  void plot(std::vector<double> &values, int numRows = 7, int numCols = 20, std::string symbol = "*");
 
   // process_mem_usage(double &, double &) - takes two doubles by reference,
   // attempts to read the system-dependent data for a process' virtual memory
   // size and resident set size, and return the results in KB.
   //
   // On failure, returns 0.0, 0.0
-  static void process_mem_usage(double& vm_usage, double& resident_set);
+  void process_mem_usage(double& vm_usage, double& resident_set);
   
-  static long getMicroTimestamp() {
+  inline long getMicroTimestamp() {
     using namespace boost::posix_time;
     ptime t0(min_date_time);
     ptime now = microsec_clock::local_time();
@@ -286,7 +303,7 @@ public:
     return td.total_microseconds();
   };
   
-  static long getRandomSeed(int param=0) {
+  inline long getRandomSeed(int param = 0) {
     using namespace boost::posix_time;
     ptime t0(min_date_time);
     ptime now = microsec_clock::local_time();
@@ -295,96 +312,97 @@ public:
     long value1 = clock();
     long value2 = random_seed_initial_value++;
 
-    return value0+value1+value2+param;
+    return value0 + value1 + value2 + param;
   };
-};