diff --git a/extensions/GeometryTools.cc b/extensions/GeometryTools.cc index f59ea8d37299506c992b909a0c2630ae0ccc1ade..45de78f22ef54ab3ab402433a3e80203cc1869ec 100644 --- a/extensions/GeometryTools.cc +++ b/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 +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 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; +} + } diff --git a/extensions/GeometryTools.h b/extensions/GeometryTools.h index 874a19b925a585f55ab06be4e3b128bb70a57430..e19b4280b6bcc940c7c76a6fd48afbcba8ef56a6 100644 --- a/extensions/GeometryTools.h +++ b/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 > &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 -void coordinate_transform(double x[], double shift[], double scale[], double alpha=0.0, double beta=0.0) -{ - mtl::dense2D 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 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 diff --git a/extensions/Helpers.h b/extensions/Helpers.h index f65b55eff1eb6fa1735d8156f24dc17119ea8a21..c5e817498ece589e018b3312331090a1728b15e9 100644 --- a/extensions/Helpers.h +++ b/extensions/Helpers.h @@ -22,44 +22,48 @@ #include #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 - 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); - }; + } template - static std::string toString(const WorldVector &value, bool brackets= true, ios_base::fmtflags flag= ios_base::scientific) + inline std::string toString(const WorldVector &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 - static std::string toString(const std::vector &vec, bool brackets= true, ios_base::fmtflags flag= ios_base::scientific) + inline std::string toString(const std::vector &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 - 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 scale, WorldVector shift, WorldVector rotate); + void transformMesh(Mesh *mesh, WorldVector scale, WorldVector shift, WorldVector rotate); - static void scaleMesh(std::vector meshes, WorldVector scale); + void scaleMesh(std::vector meshes, WorldVector scale); // scale and shift by different values in all directions - static void scaleMesh(Mesh *mesh, WorldVector shift, WorldVector scale); + void scaleMesh(Mesh *mesh, WorldVector shift, WorldVector scale); // scale by different values in all directions - static void scaleMesh(Mesh *mesh, WorldVector scale); + void scaleMesh(Mesh *mesh, WorldVector 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 getMeshDimension(Mesh *mesh); - static void getMeshDimension(Mesh *mesh, WorldVector &min_corner, WorldVector &max_corner); + WorldVector getMeshDimension(Mesh *mesh); + void getMeshDimension(Mesh *mesh, WorldVector &min_corner, WorldVector &max_corner); /// read DOFVector from AMDiS .dat-files - static void read_dof_vector(const std::string file, DOFVector *dofvec, long size); - - /// copy DOFVectors that live on different meshes - static void copyDOF(std::vector*> dof1, std::vector*> dof2, WorldVector shift, double fillValue = 0.0); - - static void copyDOF(DOFVector* dof1, DOFVector* dof2) { - WorldVector shift; shift.set(0.0); - std::vector*> dof1Vec; dof1Vec.push_back(dof1); - std::vector*> dof2Vec; dof2Vec.push_back(dof2); - copyDOF(dof1Vec, dof2Vec, shift); - }; + void read_dof_vector(const std::string file, DOFVector *dofvec, long size); /// some linear algebra methods /// =========================== - static double determinant(mtl::dense2D &A) // only for 3x3 - matrices + inline double determinant(mtl::dense2D &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 &A) // only for 3x3 - matrices + inline double determinant(WorldMatrix &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 - 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 - 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 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 - static void interpolOverLine(const Container &vec, + inline void interpolOverLine(const Container &vec, double x1, double y1, double x2, double y2, int nPoints, - std::vector > &x, std::vector &y); + std::vector > &x, std::vector &y) + { FUNCNAME("Helpers::interpolOverLine()"); + + WorldVector 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(i)/static_cast(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 *rho, std::vector, double> > &maxima); @@ -269,16 +286,16 @@ public: /// misc routines /// ============= - static void plot(std::vector &values, int numRows=7, int numCols=20, std::string symbol="*"); + void plot(std::vector &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; }; -}; - - - -/** \ingroup DOFAdministration - * \brief - * Implements a DOFIterator for a DOFIndexed object - */ -template -class DOFConstIterator : public DOFIteratorBase -{ -public: - /// Constructs a DOFIterator for cont of type t - DOFConstIterator(const DOFIndexed *obj, DOFIteratorType t) - : DOFIteratorBase(dynamic_cast(obj->getFeSpace()->getAdmin()), t), - iteratedObject(obj) - {} - - /// Constructs a DOFIterator for cont of type t - DOFConstIterator(DOFAdmin *admin, - const DOFIndexed *obj, - DOFIteratorType t) - : DOFIteratorBase(admin, t), - iteratedObject(obj) - {} - - /// Dereference operator - inline const T& operator*() - { - return *it; - } - /// Dereference operator - inline const T* operator->() - { - return &(*it); - } + +} // end namespace Helpers - inline bool operator!=(const DOFIterator& rhs) +namespace AMDiS { + + /** \ingroup DOFAdministration + * \brief + * Implements a DOFIterator for a DOFIndexed object + */ + template + class DOFConstIterator : public DOFIteratorBase { - if (this->iteratedObject != rhs.iteratedObject) - return true; + public: + /// Constructs a DOFIterator for cont of type t + DOFConstIterator(const DOFIndexed *obj, DOFIteratorType t) + : DOFIteratorBase(dynamic_cast(obj->getFeSpace()->getAdmin()), t), + iteratedObject(obj) + {} + + /// Constructs a DOFIterator for cont of type t + DOFConstIterator(DOFAdmin *admin, + const DOFIndexed *obj, + DOFIteratorType t) + : DOFIteratorBase(admin, t), + iteratedObject(obj) + {} + + /// Dereference operator + inline const T& operator*() + { + return *it; + } - if (this->it != rhs.it) - return true; + /// Dereference operator + inline const T* operator->() + { + return &(*it); + } - return false; - } + inline bool operator!=(const DOFIterator& rhs) + { + if (this->iteratedObject != rhs.iteratedObject) + return true; - inline bool operator==(const DOFIterator& rhs) - { - return !(this->operator==(rhs)); - } + if (this->it != rhs.it) + return true; -protected: - /// Implementation of DOFIteratorBase::goToBeginOfIteratedObject() - inline void goToBeginOfIteratedObject() - { - it = iteratedObject->begin(); - } + return false; + } - /// Implementation of DOFIteratorBase::goToEndOfIteratedObject() - inline void goToEndOfIteratedObject() - { - it = iteratedObject->end(); - } + inline bool operator==(const DOFIterator& rhs) + { + return !(this->operator==(rhs)); + } - /// Implementation of DOFIteratorBase::incObjectIterator() - inline void incObjectIterator() - { - ++it; - } + protected: + /// Implementation of DOFIteratorBase::goToBeginOfIteratedObject() + inline void goToBeginOfIteratedObject() + { + it = iteratedObject->begin(); + } - /// Implementation of DOFIteratorBase::incObjectIterator() - inline void decObjectIterator() - { - --it; - } + /// Implementation of DOFIteratorBase::goToEndOfIteratedObject() + inline void goToEndOfIteratedObject() + { + it = iteratedObject->end(); + } -protected: - /// Object that is iterated - const DOFIndexed *iteratedObject; + /// Implementation of DOFIteratorBase::incObjectIterator() + inline void incObjectIterator() + { + ++it; + } - /// Iterator for \ref iteratedObject - typename std::vector::const_iterator it; -}; + /// Implementation of DOFIteratorBase::incObjectIterator() + inline void decObjectIterator() + { + --it; + } + + protected: + /// Object that is iterated + const DOFIndexed *iteratedObject; -#include "Helpers.hh" + /// Iterator for \ref iteratedObject + typename std::vector::const_iterator it; + }; +} // end namespace AMDiS #endif diff --git a/extensions/Helpers.hh b/extensions/Helpers.hh deleted file mode 100644 index 2eb2377504cdd9009a08a41fad780d9ff2e3d5fa..0000000000000000000000000000000000000000 --- a/extensions/Helpers.hh +++ /dev/null @@ -1,26 +0,0 @@ - -#include "Views.h" - -template -void Helpers::interpolOverLine(const Container &vec, - double x1, double y1, double x2, double y2, - int nPoints, - std::vector > &x, std::vector &y) -{ FUNCNAME("Helpers::interpolOverLine()"); - - WorldVector p; - bool inside; - - if (nPoints <= 1) - throw(std::runtime_error("Zu wenig Punkte fuer eine Interpolation!")); - - for (int i = 0; i < nPoints; ++i) { - double lambda = static_cast(i)/static_cast(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); - } -}; diff --git a/extensions/kdtree_nanoflann_mesh.h b/extensions/kdtree_nanoflann_mesh.h index a7b7b3e0c574a6c2dc28cf86640032d1d15e3d3c..496b8e9057af1c6b1c6fdbc683654641c2d6fc53 100644 --- a/extensions/kdtree_nanoflann_mesh.h +++ b/extensions/kdtree_nanoflann_mesh.h @@ -218,14 +218,14 @@ typedef std::vector VectorOfDataType; if (getNearestElInfo(x, elInfo)) { int dim = vec.getFeSpace()->getMesh()->getDim(); DimVec lambda(dim, NO_INIT); + + double area = meshconv::triangle_area_3d(elInfo->getCoord(0).begin(), + elInfo->getCoord(1).begin(), + elInfo->getCoord(2).begin()) for (int i = 0; i < dim+1; i++) - lambda[i] = 1.0/std::max(1.e-8, norm(x - elInfo->getCoord(i))); - // normalize barycentric coords - double sum = 0.0; - for (int i = 0; i < dim+1; i++) - sum += lambda[i]; - for (int i = 0; i < dim+1; i++) - lambda[i] *= 1.0/sum; + lambda[i] = meshconv::triangle_area_3d(p.begin(), + elInfo->getCoord((i+1)%3).begin(), + elInfo->getCoord((i+2)%3).begin())/area; ElementFunctionDOFVec elFct(&vec); elFct.setElInfo(elInfo);