diff --git a/extensions/GeometryTools.cc b/extensions/GeometryTools.cc index 45de78f22ef54ab3ab402433a3e80203cc1869ec..f5338fc28d8e80d02bf0c2366608811002da82c7 100644 --- a/extensions/GeometryTools.cc +++ b/extensions/GeometryTools.cc @@ -19,6 +19,13 @@ double solve_determinant4(double x11, double x12, double x13, double x14, -x14*(x21*x32*x43 + x22*x33*x41 + x23*x31*x42 - x23*x32*x41 - x22*x31*x43 - x21*x33*x42); } +double solve_determinant3(double x11, double x12, double x13, + double x21, double x22, double x23, + double x31, double x32, double x33){ + return (x11*x22*x33 + x12*x23*x31 + x13*x21*x32) + -(x13*x22*x31 + x12*x21*x33 + x11*x23*x32); +} + // ################# // ## 2D Worlds ## // ################# @@ -1473,13 +1480,13 @@ double distance_point_triangle_3d(double point[], double tri0[], double tri1[], } -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]; +double volume_tetrahedron(double tet0[], double tet1[], double tet2[], double tet3[]){ + double det = solve_determinant3( + tet1[0]-tet0[0], tet2[0]-tet0[0], tet3[0]-tet0[0], + tet1[1]-tet0[1], tet2[1]-tet0[1], tet3[1]-tet0[1], + tet1[2]-tet0[2], tet2[2]-tet0[2], tet3[2]-tet0[2]); - return abs(Helpers::determinant(A))/6.0; + return abs(det)/6.0; } } diff --git a/extensions/GeometryTools.h b/extensions/GeometryTools.h index e19b4280b6bcc940c7c76a6fd48afbcba8ef56a6..4bee6cfef001831e66538e3ed36c41bb93ac0f6d 100644 --- a/extensions/GeometryTools.h +++ b/extensions/GeometryTools.h @@ -3,217 +3,230 @@ #ifndef GEOMETRY_TOOLBOX_H #define GEOMETRY_TOOLBOX_H -#include "AMDiS.h" +#include "FixVec.h" namespace meshconv { -// ############### -// ## General ## -// ############### +/** + * \defgroup General Dimension independet methods + * @{ + */ -//----------------------------------------< 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); - -//----------------------------------------< solve_determinant4 > +///----------------------------------------< solve_determinant4 > double solve_determinant4(double x11, double x12, double x13, double x14, double x21, double x22, double x23, double x24, double x31, double x32, double x33, double x34, double x41, double x42, double x43, double x44); -// ################# -// ## 2D Worlds ## -// ################# +///----------------------------------------< solve_determinant3 > +double solve_determinant3(double x11, double x12, double x13, + double x21, double x22, double x23, + double x31, double x32, double x33); + +/**@}*/ + +/** + * \defgroup 2d-world Method with points, lines, triangles, boxes in 2d world + * @{ + */ -//----------------------------------------< edge_length_2d > +///----------------------------------------< edge_length_2d > double edge_length_2d(double lin0[], double lin1[]); -//----------------------------------------< boundingbox_triangle_2d > +///----------------------------------------< boundingbox_triangle_2d > void boundingbox_triangle_2d(double tri0[], double tri1[], double tri2[], double min_corner[], double max_corner[]); -//----------------------------------------< centroid_of_box_2d > +///----------------------------------------< centroid_of_box_2d > void centroid_of_box_2d(double min_corner[], double max_corner[], double c[]); -//----------------------------------------< centroid_of_triangle_2d > +///----------------------------------------< centroid_of_triangle_2d > void centroid_of_triangle_2d(double tri0[], double tri1[], double tri2[], double c[]); -//----------------------------------------< point_in_triangle_2d > +///----------------------------------------< point_in_triangle_2d > bool point_in_triangle_2d(double p[], double a[], double b[], double c[]); -//----------------------------------------< point_in_triangle_generous_2d > +///----------------------------------------< point_in_triangle_generous_2d > bool point_in_triangle_generous_2d(double p[], double a[], double b[], double c[]); -//----------------------------------------< point_in_box_2d > +///----------------------------------------< point_in_box_2d > bool point_in_box_2d(double p[], double min_corner[], double max_corner[]); -//----------------------------------------< point_in_box_generous_2d > +///----------------------------------------< point_in_box_generous_2d > bool point_in_box_generous_2d(double p[], double min_corner[], double max_corner[]); -//----------------------------------------< intersection_line_line_2d > -//tests whether the line segment defined by end points 'lin0' and 'lin1' intersects the connection -//between lamp and new_dof. If the line segment is intersected in one of its vertices the index of that -//vertex (0 or 1) is returned in hit_vertex otherwise hit_vertex is set to -1. +/** + * ----------------------------------------< intersection_line_line_2d > + * tests whether the line segment defined by end points 'lin0' and 'lin1' intersects the connection + * between lamp and new_dof. If the line segment is intersected in one of its vertices the index of that + * vertex (0 or 1) is returned in hit_vertex otherwise hit_vertex is set to -1. + **/ bool intersection_line_line_2d(double lin0[], double lin1[], double lamp[], double new_dof[], long &hit_vertex); -//----------------------------------------< intersection_line_triangle_2d > +///----------------------------------------< intersection_line_triangle_2d > bool intersection_line_triangle_2d(double lin0[], double lin1[], double tri0[], double tri1[], double tri2[]); -//----------------------------------------< intersection_line_box_2d > +///----------------------------------------< intersection_line_box_2d > bool intersection_line_box_2d(double lin0[], double lin1[], double min_corner[], double max_corner[]); -//----------------------------------------< intersection_box_box_2d > +///----------------------------------------< intersection_box_box_2d > bool intersection_box_box_2d(double min_corner0[], double max_corner0[], double min_corner1[], double max_corner1[]); -//----------------------------------------< intersection_triangle_box_2d > +///----------------------------------------< intersection_triangle_box_2d > bool intersection_triangle_box_2d(double tri0[], double tri1[], double tri2[], double min_corner[], double max_corner[]); -//----------------------------------------< distance_point_line_2d > -//calculates the distance between a given point and the line segment between 'lin0' and 'lin1'. +///----------------------------------------< distance_point_line_2d > +/// calculates the distance between a given point and the line segment between 'lin0' and 'lin1'. double distance_point_line_2d(double point[], double lin0[], double lin1[]); -//----------------------------------------< distance_point_line_with_intersection_2d > -//calculates the distance between a given point and the line-segment between 'lin0' and 'lin1'. -//This version also returns the intersection-point of the perpendicular with the given line-segment. +/** + * ----------------------------------------< distance_point_line_with_intersection_2d > + * calculates the distance between a given point and the line-segment between 'lin0' and 'lin1'. + * This version also returns the intersection-point of the perpendicular with the given line-segment. + **/ double distance_point_line_with_intersection_2d(double point[], double lin0[], double lin1[], - double intersection[]); + double intersection[]); -//----------------------------------------< distance_point_box_2d > +///----------------------------------------< distance_point_box_2d > double distance_point_box_2d(double point[], double min_corner[], double max_corner[]); -//----------------------------------------< distance_point_triangle_2d > +///----------------------------------------< distance_point_triangle_2d > double distance_point_triangle_2d(double point[], double tri0[], double tri1[], double tri2[]); +///----------------------------------------< point_in_polygon > bool point_in_polygon(double point[], const std::vector<AMDiS::WorldVector<double> > &vertices); -//----------------------------------------< triangle_area_2d > +///----------------------------------------< triangle_area_2d > double triangle_area_2d(double tri0[], double tri1[], double tri2[]); -// ################# -// ## 3D Worlds ## -// ################# +/**@}*/ -//----------------------------------------< points_identical_3d > +/** + * \defgroup 3d-world Method with points, lines, triangles, boxes, terahedron in 3d world + * @{ + */ + +///----------------------------------------< points_identical_3d > bool points_identical_3d(double p1[], double p2[]); -//----------------------------------------< normal_vector_3d > +///----------------------------------------< normal_vector_3d > void normal_vector_3d(double tri0[], double tri1[], double tri2[], double normal[]); -//----------------------------------------< triangle_area_3d > +///----------------------------------------< triangle_area_3d > double triangle_area_3d(double tri0[], double tri1[], double tri2[]); -//----------------------------------------< edge_length_3d > +///----------------------------------------< edge_length_3d > double edge_length_3d(double lin0[], double lin1[]); -//----------------------------------------< triangle_max_edge_length_3d > +///----------------------------------------< triangle_max_edge_length_3d > double triangle_max_edge_length_3d(double tri0[], double tri1[], double tri2[]); -//----------------------------------------< unit_normal_vector_3d > +///----------------------------------------< unit_normal_vector_3d > void unit_normal_vector_3d(double tri0[], double tri1[], double tri2[], double normal[]); -//----------------------------------------< degenerate_triangle_3d > +///----------------------------------------< degenerate_triangle_3d > bool degenerate_triangle_3d(double tri0[], double tri1[], double tri2[]); - -//----------------------------------------< boundingbox_triangle_3d > +///----------------------------------------< boundingbox_triangle_3d > void boundingbox_triangle_3d(double tri0[], double tri1[], double tri2[], double min_corner[], double max_corner[]); -//----------------------------------------< centroid_of_box_3d > +///----------------------------------------< centroid_of_box_3d > void centroid_of_box_3d(double min_corner[], double max_corner[], double c[]); -//----------------------------------------< centroid_of_triangle_3d > +///----------------------------------------< centroid_of_triangle_3d > void centroid_of_triangle_3d(double tri0[], double tri1[], double tri2[], double c[]); -//----------------------------------------< point_in_box_3d > +///----------------------------------------< point_in_box_3d > bool point_in_box_3d(double p[], double min_corner[], double max_corner[]); -//----------------------------------------< point_in_box_generous_3d > +///----------------------------------------< point_in_box_generous_3d > bool point_in_box_generous_3d(double p[], double min_corner[], double max_corner[]); -//----------------------------------------< point_in_tetrahedron_3d > +///----------------------------------------< point_in_tetrahedron_3d > bool point_in_tetrahedron_3d(double p[], double a[], double b[], double c[], double d[]); -//----------------------------------------< intersection_line_line_3d > +///----------------------------------------< intersection_line_line_3d > bool intersection_line_line_3d(double p1[], double p2[], double p3[], double p4[]); -//----------------------------------------< intersection_line_line_with_intersection_3d > +///----------------------------------------< intersection_line_line_with_intersection_3d > bool intersection_line_line_with_intersection_3d(double p1[], double p2[], double p3[], double p4[], double intersection[]); -//----------------------------------------< intersection_line_triangle_3d_chirkov > -//ATTENTION: function does not work for lines almost in the plane of the triangle! -//calculates the intersection of a line segment and a triangle if it exists. The intersection-point is not -//calculated, the last argument "sol" remains unchanged. -// input: a line segment (lin0, lin1), and a triangle (tri0, tri1, tri2) -// return: 0 disjoint -// 1 intersection inside triangle -// 2,3,4 intersection in edge of triangle (edge_index is 0, 1, 2 respectively) -//this algorithm is based on an implementation by Nick Chirkov (journal of graphics tools 10(3):13-18, 2005). +/** + * ----------------------------------------< intersection_line_triangle_3d_chirkov > + * ATTENTION: function does not work for lines almost in the plane of the triangle! + * calculates the intersection of a line segment and a triangle if it exists. The intersection-point is not + * calculated, the last argument "sol" remains unchanged. + * input: a line segment (lin0, lin1), and a triangle (tri0, tri1, tri2) + * return: 0 disjoint + * 1 intersection inside triangle + * 2,3,4 intersection in edge of triangle (edge_index is 0, 1, 2 respectively) + * this algorithm is based on an implementation by Nick Chirkov (journal of graphics tools 10(3):13-18, 2005). + **/ int intersection_line_triangle_3d(double tri0[], double tri1[], double tri2[], - double lin0[], double lin1[], double sol[]); - -//----------------------------------------< global variables > -// int (*intersection_line_triangle_3d)(double tri0[], double tri1[], double tri2[], -// double lin0[], double lin1[], double sol[]) = intersection_line_triangle_3d_chirkov; - -//----------------------------------------< intersection_line_triangle_3d_softsurfer > -// Copyright 2001, softSurfer (www.softsurfer.com) -// This code may be freely used and modified for any purpose -// providing that this copyright notice is included with it. -// SoftSurfer makes no warranty for this code, and cannot be held -// liable for any real or imagined damage resulting from its use. -// Users of this code must verify correctness for their application. - -//calculates the intersection point of a line segment and a triangle if it exists -// input: a line segment (lin0, lin1), and a triangle (tri0, tri1, tri2) -// output: sol = intersection point (if it exists) -// return: 0 disjoint -// 1 intersection inside triangle -// 2,3,4 intersection in edge of triangle (edge_index is 0, 1, 2 respectively) + double lin0[], double lin1[], double sol[]); + +/** + * ----------------------------------------< intersection_line_triangle_3d_softsurfer > + * Copyright 2001, softSurfer (www.softsurfer.com) + * This code may be freely used and modified for any purpose + * providing that this copyright notice is included with it. + * SoftSurfer makes no warranty for this code, and cannot be held + * liable for any real or imagined damage resulting from its use. + * Users of this code must verify correctness for their application. + * + * calculates the intersection point of a line segment and a triangle if it exists + * input: a line segment (lin0, lin1), and a triangle (tri0, tri1, tri2) + * output: sol = intersection point (if it exists) + * return: 0 disjoint + * 1 intersection inside triangle + * 2,3,4 intersection in edge of triangle (edge_index is 0, 1, 2 respectively) + **/ int intersection_line_triangle_3d_softsurfer(double tri0[], double tri1[], double tri2[], - double lin0[], double lin1[], double sol[]); + double lin0[], double lin1[], double sol[]); -//----------------------------------------< intersection_line_tetrahedron_3d > +///----------------------------------------< intersection_line_tetrahedron_3d > bool intersection_line_tetrahedron_3d(double lin0[], double lin1[], double tet0[], double tet1[], double tet2[], double tet3[]); -//----------------------------------------< intersection_triangle_tetrahedron_3d > +///----------------------------------------< intersection_triangle_tetrahedron_3d > bool intersection_triangle_tetrahedron_3d(double tri0[], double tri1[], double tri2[], double tet0[], double tet1[], double tet2[], double tet3[]); -//----------------------------------------< intersection_box_box_3d > +///----------------------------------------< intersection_box_box_3d > bool intersection_box_box_3d(double min_corner0[], double max_corner0[], double min_corner1[], double max_corner1[]); -//----------------------------------------< intersection_triangle_box_3d > +///----------------------------------------< intersection_triangle_box_3d > bool intersection_triangle_box_3d(double tri0[], double tri1[], double tri2[], double min_corner[], double max_corner[]); -//----------------------------------------< distance_point_line_3d > -//calculates the distance between a given point and the line segment between 'lin0' and 'lin1'. +///----------------------------------------< distance_point_line_3d > +/// calculates the distance between a given point and the line segment between 'lin0' and 'lin1'. double distance_point_line_3d(double point[], double lin0[], double lin1[]); -//----------------------------------------< distance_point_box_3d > +///----------------------------------------< distance_point_box_3d > double distance_point_box_3d(double point[], double min_corner[], double max_corner[]); -//----------------------------------------< distance_point_triangle_3d > -//calculates the distance between a given 'point' and the triangle given by 'tri0', 'tri1', 'tri2'. +///----------------------------------------< distance_point_triangle_3d > +/// 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[]); +///----------------------------------------< coordinate_transform > template<int dow> void coordinate_transform(double x[], double shift[], double scale[], double alpha = 0.0, double beta = 0.0); -double volume_tetraeder(double tet0[], double tet1[], double tet2[], double tet3[]); +///----------------------------------------< volume_tetrahedron > +double volume_tetrahedron(double tet0[], double tet1[], double tet2[], double tet3[]); +/**@}*/ } // end namespace meshconv #include "GeometryTools.hh" diff --git a/extensions/GeometryTools.hh b/extensions/GeometryTools.hh new file mode 100644 index 0000000000000000000000000000000000000000..bfb59e01d21b3d89718f0f0c7ffe54f20ba6867f --- /dev/null +++ b/extensions/GeometryTools.hh @@ -0,0 +1,41 @@ +namespace meshconv { + + 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]; + + coords = T*coords; + + for (int i = 0; i < dow; i++) + x[i] = coords[i]+shift[i]; + } + +} // end namespace meshconv \ No newline at end of file