|
ViennaGrid - The Vienna Grid Library
2.1.0
|
00001 #ifndef VIENNAGRID_ALGORITHM_CIRCUMCENTER_HPP 00002 #define VIENNAGRID_ALGORITHM_CIRCUMCENTER_HPP 00003 00004 /* ======================================================================= 00005 Copyright (c) 2011-2014, Institute for Microelectronics, 00006 Institute for Analysis and Scientific Computing, 00007 TU Wien. 00008 00009 ----------------- 00010 ViennaGrid - The Vienna Grid Library 00011 ----------------- 00012 00013 License: MIT (X11), see file LICENSE in the base directory 00014 ======================================================================= */ 00015 00016 #include <iostream> 00017 #include <sstream> 00018 #include <string> 00019 #include <stdexcept> 00020 00021 #include "viennagrid/forwards.hpp" 00022 #include "viennagrid/algorithm/spanned_volume.hpp" 00023 #include "viennagrid/algorithm/cross_prod.hpp" 00024 #include "viennagrid/topology/all.hpp" 00025 00030 namespace viennagrid 00031 { 00032 namespace detail 00033 { 00035 template <typename ElementT, typename GeometricContainerT, typename ElementTag, typename DimensionTag> 00036 typename viennagrid::result_of::point<GeometricContainerT>::type 00037 circumcenter(ElementT const &, GeometricContainerT const &, ElementTag const &, DimensionTag const &) 00038 { 00039 typename ElementT::ERROR_COMPUTATION_OF_CIRCUMCENTER_NOT_IMPLEMENTED error_object; 00040 (void)error_object; 00041 return typename viennagrid::result_of::point<GeometricContainerT>::type(); 00042 } 00043 00044 00045 //a point is degenerate and returns its location 00047 template <typename PointAccessorT, typename ElementT> 00048 typename PointAccessorT::value_type 00049 circumcenter(PointAccessorT const accessor, ElementT const & cell, viennagrid::vertex_tag) 00050 { 00051 return accessor(cell); 00052 } 00053 00054 // 00055 // Calculation of circumcenter for a line 00056 // 00058 template <typename PointAccessorT, typename ElementT, typename DimensionTag> 00059 typename PointAccessorT::value_type 00060 circumcenter(PointAccessorT const accessor, ElementT const & cell, viennagrid::simplex_tag<1>, DimensionTag) 00061 { 00062 typedef typename PointAccessorT::value_type PointType; 00063 00064 PointType const & A = accessor( vertices(cell)[0] ); 00065 PointType const & B = accessor( vertices(cell)[1] ); 00066 00067 PointType ret = A + B; 00068 ret /= 2.0; 00069 00070 return ret; 00071 } 00072 00074 template <typename PointAccessorT, typename ElementT, typename DimensionTag> 00075 typename PointAccessorT::value_type 00076 circumcenter(PointAccessorT const accessor, ElementT const & cell, viennagrid::hypercube_tag<1>, DimensionTag) 00077 { 00078 return circumcenter(accessor, cell, viennagrid::simplex_tag<1>(), DimensionTag()); 00079 } 00080 00081 // 00082 // Calculation of circumcenter taken from Wikipedia (better reference required...) 00083 // 00085 template <typename PointAccessorT, typename ElementT> 00086 typename PointAccessorT::value_type 00087 circumcenter(PointAccessorT const accessor, ElementT const & cell, viennagrid::triangle_tag, viennagrid::dimension_tag<2>) 00088 { 00089 typedef typename PointAccessorT::value_type PointType; 00090 00091 PointType const & A = accessor( vertices(cell)[0] ); 00092 PointType const & B = accessor( vertices(cell)[1] ); 00093 PointType const & C = accessor( vertices(cell)[2] ); 00094 00095 PointType circ_cent; 00096 00097 double D = 2.0 * ( A[0] * (B[1] - C[1]) 00098 + B[0] * (C[1] - A[1]) 00099 + C[0] * (A[1] - B[1]) ); 00100 00101 double x = ( (A[1] * A[1] + A[0] * A[0]) * (B[1] - C[1]) 00102 + (B[1] * B[1] + B[0] * B[0]) * (C[1] - A[1]) 00103 + (C[1] * C[1] + C[0] * C[0]) * (A[1] - B[1]) ) / D; 00104 00105 double y = ( (A[1] * A[1] + A[0] * A[0]) * (C[0] - B[0]) 00106 + (B[1] * B[1] + B[0] * B[0]) * (A[0] - C[0]) 00107 + (C[1] * C[1] + C[0] * C[0]) * (B[0] - A[0]) ) / D; 00108 00109 circ_cent[0] = x; 00110 circ_cent[1] = y; 00111 00112 return circ_cent; 00113 } 00114 00115 00116 // 00117 // TODO: This works for rectangles only, but not for general quadrilaterals 00118 // 00120 template <typename PointAccessorT, typename ElementT> 00121 typename PointAccessorT::value_type 00122 circumcenter(PointAccessorT const accessor, ElementT const & cell, viennagrid::quadrilateral_tag, viennagrid::dimension_tag<2>) 00123 { 00124 typedef typename PointAccessorT::value_type PointType; 00125 typedef typename viennagrid::result_of::coord<PointType>::type CoordType; 00126 00127 typedef typename viennagrid::result_of::const_element_range<ElementT, vertex_tag>::type VertexOnCellRange; 00128 typedef typename viennagrid::result_of::iterator<VertexOnCellRange>::type VertexOnCellIterator; 00129 00130 PointType p0(0.0, 0.0); 00131 00132 VertexOnCellRange vertices = viennagrid::elements<vertex_tag>(cell); 00133 for (VertexOnCellIterator vocit = vertices.begin(); 00134 vocit != vertices.end(); 00135 ++vocit) 00136 { 00137 p0 += accessor( *vocit ); 00138 } 00139 00140 p0 /= static_cast<CoordType>(vertices.size()); 00141 00142 return p0; 00143 } 00144 00145 00146 00147 // 00148 // Calculation of circumcenter taken from Wikipedia (better reference required...) 00149 // 00151 template <typename PointAccessorT, typename ElementT> 00152 typename PointAccessorT::value_type 00153 circumcenter(PointAccessorT const accessor, ElementT const & cell, viennagrid::triangle_tag, viennagrid::dimension_tag<3>) 00154 { 00155 typedef typename PointAccessorT::value_type PointType; 00156 00157 PointType const & A = accessor( vertices(cell)[0] ); 00158 PointType const & B = accessor( vertices(cell)[1] ); 00159 PointType const & C = accessor( vertices(cell)[2] ); 00160 00161 double denominator = 2.0 * viennagrid::inner_prod(viennagrid::cross_prod(A-B, B-C), viennagrid::cross_prod(A-B, B-C)); 00162 00163 double alpha = viennagrid::inner_prod(B - C, B - C) * viennagrid::inner_prod(A - B, A - C); 00164 double beta = viennagrid::inner_prod(A - C, A - C) * viennagrid::inner_prod(B - A, B - C); 00165 double gamma = viennagrid::inner_prod(A - B, A - B) * viennagrid::inner_prod(C - A, C - B); 00166 00167 00168 PointType circ_cent = A * (alpha/denominator) + B * (beta/denominator) + C * (gamma/denominator); 00169 00170 return circ_cent; 00171 } 00172 00173 00174 // 00175 // Note: This works for rectangles only, but not for general quadrilaterals 00176 // 00178 template <typename PointAccessorT, typename ElementT> 00179 typename PointAccessorT::value_type 00180 circumcenter(PointAccessorT const accessor, ElementT const & cell, viennagrid::quadrilateral_tag, viennagrid::dimension_tag<3>) 00181 { 00182 typedef typename PointAccessorT::value_type PointType; 00183 typedef typename viennagrid::result_of::coord<PointType>::type CoordType; 00184 00185 typedef typename viennagrid::result_of::const_element_range<ElementT, vertex_tag>::type VertexOnCellRange; 00186 typedef typename viennagrid::result_of::iterator<VertexOnCellRange>::type VertexOnCellIterator; 00187 00188 PointType p0(0.0, 0.0, 0.0); 00189 00190 VertexOnCellRange vertices = viennagrid::elements<vertex_tag>(cell); 00191 for (VertexOnCellIterator vocit = vertices.begin(); 00192 vocit != vertices.end(); 00193 ++vocit) 00194 { 00195 p0 += accessor( *vocit ); 00196 } 00197 00198 p0 /= static_cast<CoordType>(vertices.size()); 00199 00200 return p0; 00201 } 00202 00203 00204 // 00205 // Calculation of circumcenter taken from Wikipedia (better reference required...) 00206 // 00208 template <typename PointAccessorT, typename ElementT> 00209 typename PointAccessorT::value_type 00210 circumcenter(PointAccessorT const accessor, ElementT const & cell, viennagrid::tetrahedron_tag, viennagrid::dimension_tag<3>) 00211 { 00212 typedef typename PointAccessorT::value_type PointType; 00213 00214 typedef typename viennagrid::result_of::const_element_range<ElementT, vertex_tag>::type VertexOnCellRange; 00215 typedef typename viennagrid::result_of::iterator<VertexOnCellRange>::type VertexOnCellIterator; 00216 00217 VertexOnCellRange vertices = viennagrid::elements<vertex_tag>(cell); 00218 VertexOnCellIterator vocit = vertices.begin(); 00219 00220 PointType const & O = accessor(*vocit); ++vocit; 00221 PointType const & A = accessor(*vocit) - O; ++vocit; 00222 PointType const & B = accessor(*vocit) - O; ++vocit; 00223 PointType const & C = accessor(*vocit) - O; 00224 00225 PointType circ_cent = (cross_prod(A, B) * viennagrid::inner_prod(C, C) 00226 + cross_prod(C, A) * viennagrid::inner_prod(B, B) 00227 + cross_prod(B, C) * viennagrid::inner_prod(A, A) 00228 ) / (viennagrid::inner_prod(A, viennagrid::cross_prod(B, C)) * 2.0); 00229 PointType cprodBC = viennagrid::cross_prod(B, C); 00230 if (std::fabs(viennagrid::inner_prod(A, cprodBC)) < 1e-10 * viennagrid::norm(cprodBC) ) 00231 { 00232 std::cerr << "Near singularity in circum center calculation!" << std::endl; 00233 std::cerr << "A: " << A << std::endl; 00234 std::cerr << "B: " << B << std::endl; 00235 std::cerr << "C: " << C << std::endl; 00236 std::cerr << "B x C: " << viennagrid::cross_prod(B, C) << std::endl; 00237 throw std::runtime_error("Near singularity in circum center calculation!"); 00238 } 00239 return circ_cent + O; 00240 } 00241 00242 // 00243 // Note: This works for rectangles only, but not for general quadrilaterals 00244 // 00246 template <typename PointAccessorT, typename ElementT> 00247 typename PointAccessorT::value_type 00248 circumcenter(PointAccessorT const accessor, ElementT const & cell, viennagrid::hexahedron_tag, viennagrid::dimension_tag<3>) 00249 { 00250 typedef typename PointAccessorT::value_type PointType; 00251 typedef typename viennagrid::result_of::coord<PointType>::type CoordType; 00252 00253 typedef typename viennagrid::result_of::const_element_range<ElementT, vertex_tag>::type VertexOnCellRange; 00254 typedef typename viennagrid::result_of::iterator<VertexOnCellRange>::type VertexOnCellIterator; 00255 00256 PointType p0(0.0, 0.0); 00257 00258 VertexOnCellRange vertices = viennagrid::elements<vertex_tag>(cell); 00259 for (VertexOnCellIterator vocit = vertices.begin(); 00260 vocit != vertices.end(); 00261 ++vocit) 00262 { 00263 p0 += accessor(*vocit); 00264 } 00265 00266 p0 /= static_cast<CoordType>(vertices.size()); 00267 00268 return p0; 00269 } 00270 00271 } //namespace detail 00272 00278 template <typename PointAccessorT, typename ElementT> 00279 typename PointAccessorT::value_type 00280 circumcenter(PointAccessorT const accessor, ElementT const & element) 00281 { 00282 typedef typename PointAccessorT::value_type PointType; 00283 00284 return detail::circumcenter(accessor, 00285 element, 00286 typename ElementT::tag(), 00287 viennagrid::dimension_tag< result_of::static_size<PointType>::value >()); 00288 } 00289 00290 00295 template<typename ElementT> 00296 typename viennagrid::result_of::point< ElementT >::type 00297 circumcenter(ElementT const & element) 00298 { 00299 return circumcenter( default_point_accessor(element), element ); 00300 } 00301 00302 00303 } //namespace viennagrid 00304 #endif
1.7.6.1