00001 #ifndef VIENNAGRID_ALGORITHM_CIRCUMCENTER_HPP
00002 #define VIENNAGRID_ALGORITHM_CIRCUMCENTER_HPP
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021 #include <iostream>
00022 #include <sstream>
00023 #include <string>
00024 #include <stdexcept>
00025
00026 #include "viennagrid/forwards.h"
00027 #include "viennagrid/algorithm/spanned_volume.hpp"
00028 #include "viennagrid/algorithm/cross_prod.hpp"
00029
00030 #include "viennadata/api.hpp"
00031
00036 namespace viennagrid
00037 {
00038 namespace detail
00039 {
00041 template <typename ElementType, typename ElementTag, typename DimensionTag>
00042 typename viennagrid::result_of::point<typename ElementType::config_type>::type
00043 circumcenter(ElementType const & cell, ElementTag const &, DimensionTag const &)
00044 {
00045 typedef typename ElementType::ERROR_COMPUTATION_OF_CIRCUMCENTER_NOT_IMPLEMENTED error_type;
00046 return typename viennagrid::result_of::point<typename ElementType::config_type>::type();
00047 }
00048
00049
00050
00052 template <typename ElementType>
00053 typename viennagrid::result_of::point<typename ElementType::config_type>::type
00054 circumcenter(ElementType const & cell, viennagrid::point_tag)
00055 {
00056 return cell.point();
00057 }
00058
00059
00060
00061
00063 template <typename ElementType, typename DimensionTag>
00064 typename viennagrid::result_of::point<typename ElementType::config_type>::type
00065 circumcenter(ElementType const & cell, viennagrid::simplex_tag<1>, DimensionTag)
00066 {
00067 typedef typename ElementType::config_type Config;
00068 typedef typename viennagrid::result_of::point<Config>::type PointType;
00069
00070 typedef typename viennagrid::result_of::const_ncell_range<ElementType, 0>::type VertexOnCellRange;
00071
00072 VertexOnCellRange vertices = viennagrid::ncells<0>(cell);
00073
00074 PointType const & A = vertices[0].point();
00075 PointType const & B = vertices[1].point();
00076
00077 PointType ret = A + B;
00078 ret /= 2.0;
00079
00080 return ret;
00081 }
00082
00084 template <typename ElementType, typename DimensionTag>
00085 typename viennagrid::result_of::point<typename ElementType::config_type>::type
00086 circumcenter(ElementType const & cell, viennagrid::hypercube_tag<1>, DimensionTag)
00087 {
00088 return circumcenter(cell, viennagrid::simplex_tag<1>(), DimensionTag());
00089 }
00090
00091
00092
00093
00095 template <typename ElementType>
00096 typename viennagrid::result_of::point<typename ElementType::config_type>::type
00097 circumcenter(ElementType const & cell, viennagrid::triangle_tag, viennagrid::dimension_tag<2>)
00098 {
00099 typedef typename ElementType::config_type Config;
00100 typedef typename Config::cell_tag CellTag;
00101 typedef typename viennagrid::result_of::point<Config>::type PointType;
00102 typedef typename viennagrid::result_of::ncell<Config, 0>::type VertexType;
00103 typedef typename viennagrid::result_of::ncell<Config, 1>::type EdgeType;
00104
00105 typedef typename viennagrid::result_of::const_ncell_range<ElementType, 0>::type VertexOnCellRange;
00106
00107 VertexOnCellRange vertices = viennagrid::ncells<0>(cell);
00108
00109 PointType const & A = vertices[0].point();
00110 PointType const & B = vertices[1].point();
00111 PointType const & C = vertices[2].point();
00112
00113 PointType circ_cent;
00114
00115 double D = 2.0 * ( A[0] * (B[1] - C[1])
00116 + B[0] * (C[1] - A[1])
00117 + C[0] * (A[1] - B[1]) );
00118
00119 double x = ( (A[1] * A[1] + A[0] * A[0]) * (B[1] - C[1])
00120 + (B[1] * B[1] + B[0] * B[0]) * (C[1] - A[1])
00121 + (C[1] * C[1] + C[0] * C[0]) * (A[1] - B[1]) ) / D;
00122
00123 double y = ( (A[1] * A[1] + A[0] * A[0]) * (C[0] - B[0])
00124 + (B[1] * B[1] + B[0] * B[0]) * (A[0] - C[0])
00125 + (C[1] * C[1] + C[0] * C[0]) * (B[0] - A[0]) ) / D;
00126
00127 circ_cent[0] = x;
00128 circ_cent[1] = y;
00129
00130 return circ_cent;
00131 }
00132
00133
00134
00135
00136
00138 template <typename CellType>
00139 typename viennagrid::result_of::point<typename CellType::config_type>::type
00140 circumcenter(CellType const & cell, viennagrid::quadrilateral_tag, viennagrid::dimension_tag<2>)
00141 {
00142 typedef typename CellType::config_type Config;
00143 typedef typename Config::cell_tag CellTag;
00144 typedef typename viennagrid::result_of::point<Config>::type PointType;
00145 typedef typename viennagrid::result_of::ncell<Config, 0>::type VertexType;
00146 typedef typename viennagrid::result_of::ncell<Config, 1>::type EdgeType;
00147
00148 typedef typename viennagrid::result_of::const_ncell_range<CellType, 0>::type VertexOnCellRange;
00149 typedef typename viennagrid::result_of::iterator<VertexOnCellRange>::type VertexOnCellIterator;
00150
00151 PointType p0(0.0, 0.0);
00152
00153 VertexOnCellRange vertices = viennagrid::ncells<0>(cell);
00154 for (VertexOnCellIterator vocit = vertices.begin();
00155 vocit != vertices.end();
00156 ++vocit)
00157 {
00158 p0 += vocit->point();
00159 }
00160
00161 p0 /= viennagrid::topology::bndcells<CellTag, 0>::num;
00162
00163 return p0;
00164 }
00165
00166
00167
00168
00169
00170
00172 template <typename ElementType>
00173 typename viennagrid::result_of::point<typename ElementType::config_type>::type
00174 circumcenter(ElementType const & cell, viennagrid::triangle_tag, viennagrid::dimension_tag<3>)
00175 {
00176 typedef typename ElementType::config_type Config;
00177 typedef typename Config::cell_tag CellTag;
00178 typedef typename viennagrid::result_of::point<Config>::type PointType;
00179 typedef typename viennagrid::result_of::ncell<Config, 0>::type VertexType;
00180 typedef typename viennagrid::result_of::ncell<Config, 1>::type EdgeType;
00181
00182 typedef typename viennagrid::result_of::const_ncell_range<ElementType, 0>::type VertexOnCellRange;
00183
00184 VertexOnCellRange vertices = viennagrid::ncells<0>(cell);
00185
00186 PointType const & A = vertices[0].point();
00187 PointType const & B = vertices[1].point();
00188 PointType const & C = vertices[2].point();
00189
00190 double denominator = 2.0 * viennagrid::inner_prod(viennagrid::cross_prod(A-B, B-C), viennagrid::cross_prod(A-B, B-C));
00191
00192 double alpha = viennagrid::inner_prod(B - C, B - C) * viennagrid::inner_prod(A - B, A - C);
00193 double beta = viennagrid::inner_prod(A - C, A - C) * viennagrid::inner_prod(B - A, B - C);
00194 double gamma = viennagrid::inner_prod(A - B, A - B) * viennagrid::inner_prod(C - A, C - B);
00195
00196
00197 PointType circ_cent = A * (alpha/denominator) + B * (beta/denominator) + C * (gamma/denominator);
00198
00199 return circ_cent;
00200 }
00201
00202
00203
00204
00205
00207 template <typename CellType>
00208 typename viennagrid::result_of::point<typename CellType::config_type>::type
00209 circumcenter(CellType const & cell, viennagrid::quadrilateral_tag, viennagrid::dimension_tag<3>)
00210 {
00211 typedef typename CellType::config_type Config;
00212 typedef typename CellType::tag ElementTag;
00213 typedef typename viennagrid::result_of::point<Config>::type PointType;
00214 typedef typename viennagrid::result_of::ncell<Config, 0>::type VertexType;
00215 typedef typename viennagrid::result_of::ncell<Config, 1>::type EdgeType;
00216
00217 typedef typename viennagrid::result_of::const_ncell_range<CellType, 0>::type VertexOnCellRange;
00218 typedef typename viennagrid::result_of::iterator<VertexOnCellRange>::type VertexOnCellIterator;
00219
00220 PointType p0(0.0, 0.0, 0.0);
00221
00222 VertexOnCellRange vertices = viennagrid::ncells<0>(cell);
00223 for (VertexOnCellIterator vocit = vertices.begin();
00224 vocit != vertices.end();
00225 ++vocit)
00226 {
00227 p0 += vocit->point();
00228 }
00229
00230 p0 /= viennagrid::topology::bndcells<ElementTag, 0>::num;
00231
00232 return p0;
00233 }
00234
00235
00236
00237
00238
00240 template <typename ElementType>
00241 typename viennagrid::result_of::point<typename ElementType::config_type>::type
00242 circumcenter(ElementType const & cell, viennagrid::tetrahedron_tag, viennagrid::dimension_tag<3>)
00243 {
00244 typedef typename ElementType::config_type Config;
00245 typedef typename Config::cell_tag CellTag;
00246 typedef typename viennagrid::result_of::point<Config>::type PointType;
00247 typedef typename viennagrid::result_of::ncell<Config, 0>::type VertexType;
00248 typedef typename viennagrid::result_of::ncell<Config, 1>::type EdgeType;
00249
00250 typedef typename viennagrid::result_of::const_ncell_range<ElementType, 0>::type VertexOnCellRange;
00251 typedef typename viennagrid::result_of::iterator<VertexOnCellRange>::type VertexOnCellIterator;
00252
00253 VertexOnCellRange vertices = viennagrid::ncells<0>(cell);
00254 VertexOnCellIterator vocit = vertices.begin();
00255
00256 PointType const & O = vocit->point(); ++vocit;
00257 PointType const & A = vocit->point() - O; ++vocit;
00258 PointType const & B = vocit->point() - O; ++vocit;
00259 PointType const & C = vocit->point() - O;
00260
00261 PointType circ_cent = (cross_prod(A, B) * viennagrid::inner_prod(C, C)
00262 + cross_prod(C, A) * viennagrid::inner_prod(B, B)
00263 + cross_prod(B, C) * viennagrid::inner_prod(A, A)
00264 ) / (viennagrid::inner_prod(A, viennagrid::cross_prod(B, C)) * 2.0);
00265 if (viennagrid::inner_prod(A, viennagrid::cross_prod(B, C)) == 0)
00266 {
00267 std::cout << "Singularity in circum center calculation!" << std::endl;
00268 std::cout << "A: " << A << std::endl;
00269 std::cout << "B: " << B << std::endl;
00270 std::cout << "C: " << C << std::endl;
00271 std::cout << "B x C: " << viennagrid::cross_prod(B, C) << std::endl;
00272 exit(0);
00273 }
00274 return circ_cent + O;
00275 }
00276
00277
00278
00279
00281 template <typename CellType>
00282 typename viennagrid::result_of::point<typename CellType::config_type>::type
00283 circumcenter(CellType const & cell, viennagrid::hexahedron_tag, viennagrid::dimension_tag<3>)
00284 {
00285 typedef typename CellType::config_type Config;
00286 typedef typename Config::cell_tag CellTag;
00287 typedef typename viennagrid::result_of::point<Config>::type PointType;
00288 typedef typename viennagrid::result_of::ncell<Config, 0>::type VertexType;
00289 typedef typename viennagrid::result_of::ncell<Config, 1>::type EdgeType;
00290
00291 typedef typename viennagrid::result_of::const_ncell_range<CellType, 0>::type VertexOnCellRange;
00292 typedef typename viennagrid::result_of::iterator<VertexOnCellRange>::type VertexOnCellIterator;
00293
00294 PointType p0(0.0, 0.0);
00295
00296 VertexOnCellRange vertices = viennagrid::ncells<0>(cell);
00297 for (VertexOnCellIterator vocit = vertices.begin();
00298 vocit != vertices.end();
00299 ++vocit)
00300 {
00301 p0 += vocit->point();
00302 }
00303
00304 p0 /= viennagrid::topology::bndcells<CellTag, 0>::num;
00305
00306 return p0;
00307 }
00308
00309 }
00310
00315 template <typename CellType>
00316 typename viennagrid::result_of::point<typename CellType::config_type>::type
00317 circumcenter(CellType const & cell)
00318 {
00319 return detail::circumcenter(cell,
00320 typename CellType::tag(),
00321 viennagrid::dimension_tag<CellType::config_type::coordinate_system_tag::dim>());
00322 }
00323
00324 }
00325 #endif