00001 #ifndef VIENNAGRID_ALGORITHM_SPANNED_VOLUME_HPP
00002 #define VIENNAGRID_ALGORITHM_SPANNED_VOLUME_HPP
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023 #include <math.h>
00024 #include "viennagrid/forwards.h"
00025 #include "viennagrid/algorithm/cross_prod.hpp"
00026 #include "viennagrid/algorithm/norm.hpp"
00027 #include "viennagrid/algorithm/inner_prod.hpp"
00028 #include "viennagrid/traits/point.hpp"
00029
00030
00035 namespace viennagrid
00036 {
00037
00038 namespace detail
00039 {
00040 template <typename PointType,
00041 long dim = traits::dimension<PointType>::value>
00042 struct spanned_volume_impl;
00043
00045 template <typename PointType>
00046 struct spanned_volume_impl<PointType, 1>
00047 {
00048 typedef typename traits::value_type<PointType>::type value_type;
00049
00050 static value_type apply(PointType const & p1,
00051 PointType const & p2)
00052 {
00053 return fabs(p2[0] - p1[0]);
00054 }
00055 };
00056
00057
00059 template <typename PointType>
00060 struct spanned_volume_impl<PointType, 2>
00061 {
00062 typedef typename traits::value_type<PointType>::type value_type;
00063
00064 static value_type apply(PointType const & p1,
00065 PointType const & p2)
00066 {
00067
00068 return sqrt( (p2[0] - p1[0]) * (p2[0] - p1[0])
00069 + (p2[1] - p1[1]) * (p2[1] - p1[1]) );
00070 }
00071
00072 static value_type apply(PointType const & A,
00073 PointType const & B,
00074 PointType const & C)
00075 {
00076
00077 return fabs( A[0] * (B[1] - C[1])
00078 + B[0] * (C[1] - A[1])
00079 + C[0] * (A[1] - B[1]) ) / 2.0;
00080 }
00081
00082 };
00083
00084
00086 template <typename PointType>
00087 struct spanned_volume_impl<PointType, 3>
00088 {
00089 typedef typename traits::value_type<PointType>::type value_type;
00090
00091 static value_type apply(PointType const & p1,
00092 PointType const & p2)
00093 {
00094
00095 return sqrt( (p2[0] - p1[0]) * (p2[0] - p1[0])
00096 + (p2[1] - p1[1]) * (p2[1] - p1[1])
00097 + (p2[2] - p1[2]) * (p2[2] - p1[2]) );
00098 }
00099
00100 static value_type apply(PointType const & p1,
00101 PointType const & p2,
00102 PointType const & p3)
00103 {
00104 PointType v1 = p2 - p1;
00105 PointType v2 = p3 - p1;
00106
00107 PointType v3 = cross_prod(v1, v2);
00108
00109 return norm(v3) / 2.0;
00110 }
00111
00112 static value_type apply(PointType const & p1,
00113 PointType const & p2,
00114 PointType const & p3,
00115 PointType const & p4)
00116 {
00117 PointType v1 = p2 - p1;
00118 PointType v2 = p3 - p1;
00119 PointType v3 = p4 - p1;
00120
00121 return fabs(inner_prod(v1, cross_prod(v2, v3)) ) / 6.0;
00122 }
00123
00124 };
00125 }
00126
00127
00128
00129
00130
00132 template<typename PointType1, typename PointType2, typename CSystem1, typename CSystem2>
00133 typename traits::value_type<PointType1>::type
00134 spanned_volume_impl(PointType1 const & p1,
00135 PointType2 const & p2,
00136 CSystem1 const &,
00137 CSystem2 const &)
00138 {
00139 typedef typename traits::value_type<PointType1>::type value_type;
00140 typedef typename result_of::cartesian_point<PointType1>::type CartesianPoint1;
00141
00142 return detail::spanned_volume_impl<CartesianPoint1>::apply(to_cartesian(p1), to_cartesian(p2));
00143 }
00144
00146 template<typename PointType1, typename PointType2, typename PointType3,
00147 typename CSystem1, typename CSystem2, typename CSystem3>
00148 typename traits::value_type<PointType1>::type
00149 spanned_volume_impl(PointType1 const & p1,
00150 PointType2 const & p2,
00151 PointType3 const & p3,
00152 CSystem1 const &,
00153 CSystem2 const &,
00154 CSystem3 const &)
00155 {
00156 typedef typename traits::value_type<PointType1>::type value_type;
00157 typedef typename result_of::cartesian_point<PointType1>::type CartesianPoint1;
00158
00159 return detail::spanned_volume_impl<CartesianPoint1>::apply(to_cartesian(p1), to_cartesian(p2), to_cartesian(p3));
00160 }
00161
00163 template<typename PointType1, typename PointType2, typename PointType3, typename PointType4,
00164 typename CSystem1, typename CSystem2, typename CSystem3, typename CSystem4>
00165 typename traits::value_type<PointType1>::type
00166 spanned_volume_impl(PointType1 const & p1,
00167 PointType2 const & p2,
00168 PointType3 const & p3,
00169 PointType4 const & p4,
00170 CSystem1 const &,
00171 CSystem2 const &,
00172 CSystem3 const &,
00173 CSystem4 const &)
00174 {
00175 typedef typename traits::value_type<PointType1>::type value_type;
00176 typedef typename result_of::cartesian_point<PointType1>::type CartesianPoint1;
00177
00178 return detail::spanned_volume_impl<CartesianPoint1>::apply(to_cartesian(p1), to_cartesian(p2), to_cartesian(p3), to_cartesian(p4));
00179 }
00180
00181
00182
00183
00185 template<typename PointType1, typename PointType2, long d>
00186 typename traits::value_type<PointType1>::type
00187 spanned_volume_impl(PointType1 const & p1,
00188 PointType2 const & p2,
00189 cartesian_cs<d>,
00190 cartesian_cs<d>)
00191 {
00192 return detail::spanned_volume_impl<PointType1>::apply(p1, p2);
00193 }
00194
00196 template <typename PointType1, typename PointType2, typename PointType3, long d>
00197 typename traits::value_type<PointType1>::type
00198 spanned_volume_impl(PointType1 const & p1,
00199 PointType2 const & p2,
00200 PointType3 const & p3,
00201 cartesian_cs<d>,
00202 cartesian_cs<d>,
00203 cartesian_cs<d>)
00204 {
00205 return detail::spanned_volume_impl<PointType1>::apply(p1, p2, p3);
00206 }
00207
00209 template <typename PointType1, typename PointType2, typename PointType3, typename PointType4, long d>
00210 typename traits::value_type<PointType1>::type
00211 spanned_volume_impl(PointType1 const & p1,
00212 PointType2 const & p2,
00213 PointType3 const & p3,
00214 PointType4 const & p4,
00215 cartesian_cs<d>,
00216 cartesian_cs<d>,
00217 cartesian_cs<d>,
00218 cartesian_cs<d>)
00219 {
00220 return detail::spanned_volume_impl<PointType1>::apply(p1, p2, p3, p4);
00221 }
00222
00223
00224
00225
00226
00227
00229 template <typename PointType1, typename PointType2>
00230 typename traits::value_type<PointType1>::type
00231 spanned_volume(PointType1 const & p1, PointType2 const & p2)
00232 {
00233 return spanned_volume_impl(p1,
00234 p2,
00235 typename traits::coordinate_system<PointType1>::type(),
00236 typename traits::coordinate_system<PointType2>::type());
00237 }
00238
00239
00241 template <typename PointType1, typename PointType2, typename PointType3>
00242 typename traits::value_type<PointType1>::type
00243 spanned_volume(PointType1 const & p1, PointType2 const & p2, PointType3 const & p3)
00244 {
00245 return spanned_volume_impl(p1,
00246 p2,
00247 p3,
00248 typename traits::coordinate_system<PointType1>::type(),
00249 typename traits::coordinate_system<PointType2>::type(),
00250 typename traits::coordinate_system<PointType3>::type()
00251 );
00252
00253 }
00254
00255
00257 template <typename PointType1, typename PointType2, typename PointType3, typename PointType4>
00258 typename traits::value_type<PointType1>::type
00259 spanned_volume(PointType1 const & p1,
00260 PointType2 const & p2,
00261 PointType3 const & p3,
00262 PointType4 const & p4)
00263 {
00264 return spanned_volume_impl(p1,
00265 p2,
00266 p3,
00267 p4,
00268 typename traits::coordinate_system<PointType1>::type(),
00269 typename traits::coordinate_system<PointType2>::type(),
00270 typename traits::coordinate_system<PointType3>::type(),
00271 typename traits::coordinate_system<PointType4>::type()
00272 );
00273 }
00274
00275
00276
00277 }
00278 #endif