• Main Page
  • Namespaces
  • Data Structures
  • Files
  • File List
  • Globals

/export/development/ViennaGrid/release/ViennaGrid-1.0.0/viennagrid/algorithm/circumcenter.hpp

Go to the documentation of this file.
00001 #ifndef VIENNAGRID_ALGORITHM_CIRCUMCENTER_HPP
00002 #define VIENNAGRID_ALGORITHM_CIRCUMCENTER_HPP
00003 
00004 /* =======================================================================
00005    Copyright (c) 2011, Institute for Microelectronics,
00006                        Institute for Analysis and Scientific Computing,
00007                        TU Wien.
00008 
00009                             -----------------
00010                      ViennaGrid - The Vienna Grid Library
00011                             -----------------
00012 
00013    Authors:      Karl Rupp                           rupp@iue.tuwien.ac.at
00014                  Josef Weinbub                    weinbub@iue.tuwien.ac.at
00015                
00016    (A list of additional contributors can be found in the PDF manual)
00017 
00018    License:      MIT (X11), see file LICENSE in the base directory
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     //a point is degenerate and returns its location
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     // Calculation of circumcenter for a line
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     // Calculation of circumcenter taken from Wikipedia (better reference required...)
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     // TODO: This works for rectangles only, but not for general quadrilaterals
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     // Calculation of circumcenter taken from Wikipedia (better reference required...)
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     // Note: This works for rectangles only, but not for general quadrilaterals
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     // Calculation of circumcenter taken from Wikipedia (better reference required...)
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     // Note: This works for rectangles only, but not for general quadrilaterals
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   } //namespace detail
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 } //namespace viennagrid
00325 #endif

Generated on Wed Sep 14 2011 19:21:30 for ViennaGrid by  doxygen 1.7.1