Go to the documentation of this file.00001 #ifndef VIENNAMATH_RUNTIME_NUMERICAL_QUADRATURE_HPP
00002 #define VIENNAMATH_RUNTIME_NUMERICAL_QUADRATURE_HPP
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020 #include <string>
00021 #include <vector>
00022 #include <math.h>
00023 #include <iostream>
00024 #include <sstream>
00025 #include "viennamath/forwards.h"
00026 #include "viennamath/runtime/integral.hpp"
00027 #include "viennamath/manipulation/substitute.hpp"
00028 #include "viennamath/manipulation/eval.hpp"
00029
00034 namespace viennamath
00035 {
00036
00038 template <typename InterfaceType>
00039 class numerical_quadrature_interface
00040 {
00041 public:
00042 typedef typename InterfaceType::numeric_type numeric_type;
00043
00045 virtual numerical_quadrature_interface * clone() const = 0;
00046
00048 virtual numeric_type eval(rt_interval<InterfaceType> const & interv,
00049 rt_expr<InterfaceType> const & e,
00050 rt_variable<InterfaceType> const & var) const = 0;
00051 };
00052
00054 template <typename InterfaceType = default_interface_type>
00055 class rt_gauss_quad_1 : public numerical_quadrature_interface<InterfaceType>
00056 {
00057 typedef numerical_quadrature_interface<InterfaceType> BaseType;
00058
00059 public:
00060 typedef typename BaseType::numeric_type numeric_type;
00061
00062 enum { id = 1 };
00063
00064 BaseType * clone() const { return new rt_gauss_quad_1(); }
00065
00066
00067 numeric_type eval(rt_interval<InterfaceType> const & interv,
00068 rt_expr<InterfaceType> const & e,
00069 rt_variable<InterfaceType> const & var) const
00070 {
00071 assert(interv.lower().get()->is_constant() &&
00072 interv.upper().get()->is_constant() &&
00073 "Upper or lower interval is not a constant!");
00074
00075 numeric_type abscissa = (viennamath::eval(interv.lower(), 0.0) + viennamath::eval(interv.upper(), 0.0)) / 2.0;
00076 rt_expr<InterfaceType> tmp = viennamath::substitute(var, abscissa, e);
00077 return tmp.get()->unwrap() * viennamath::eval(interv.upper() - interv.lower(), 0.0);
00078 }
00079 private:
00080 std::vector<numeric_type> p_;
00081 };
00082 typedef rt_gauss_quad_1<> gauss_quad_1;
00083
00084
00085
00086
00087
00089 template <typename InterfaceType = default_interface_type>
00090 class rt_numerical_quadrature
00091 {
00092 public:
00093 typedef typename InterfaceType::numeric_type NumericT;
00094
00095
00096
00097
00098
00099
00100
00101
00102
00103
00104
00105
00106
00107 rt_numerical_quadrature() {}
00108
00109
00110
00111
00112
00113 rt_numerical_quadrature(numerical_quadrature_interface<InterfaceType> * ptr) : quadrature_rule_(ptr) {}
00114
00115
00116
00117 rt_numerical_quadrature(const rt_numerical_quadrature & other) : quadrature_rule_(other.quadrature_rule_->clone()) {}
00118
00119
00120
00121 rt_numerical_quadrature & operator=(const rt_numerical_quadrature & other)
00122 {
00123 quadrature_rule_ = std::auto_ptr<numerical_quadrature_interface<InterfaceType> >(other.quadrature_rule_->clone());
00124 return *this;
00125 }
00126
00127
00129 NumericT operator()(rt_expr<InterfaceType> const & e) const
00130 {
00131 const rt_unary_expr<InterfaceType> * integral_expression = dynamic_cast<const rt_unary_expr<InterfaceType> *>(e.get());
00132
00133 if (integral_expression != NULL)
00134 {
00135 typedef op_unary<op_rt_integral<InterfaceType>, InterfaceType> IntegrationOperation;
00136
00137 const IntegrationOperation * op_tmp = dynamic_cast<const IntegrationOperation *>(integral_expression->op());
00138
00139 if (op_tmp != NULL)
00140 {
00141 return quadrature_rule_->eval(op_tmp->op().interval(),
00142 rt_expr<InterfaceType>(integral_expression->lhs()->clone()),
00143 op_tmp->op().variable());
00144 }
00145 else
00146 {
00147 typedef op_unary<op_rt_symbolic_integral<InterfaceType>, InterfaceType> SymbolicIntegrationOperation;
00148
00149 const SymbolicIntegrationOperation * op_symb_tmp = dynamic_cast<const SymbolicIntegrationOperation *>(integral_expression->op());
00150
00151 if (op_symb_tmp != NULL)
00152 {
00153 return quadrature_rule_->eval(op_tmp->op().interval(),
00154 rt_expr<InterfaceType>(integral_expression->lhs()->clone()),
00155 op_tmp->op().variable());
00156 }
00157
00158
00159
00160 throw integration_without_integral_exception();
00161 }
00162 }
00163 else
00164 {
00165 if (e.get()->is_constant())
00166 {
00167 if (e.get()->eval(0) == 0)
00168 return rt_constant<typename InterfaceType::numeric_type, InterfaceType>(0);
00169 }
00170
00171
00172
00173 throw integration_without_integral_exception();
00174 }
00175
00176 return 0;
00177 }
00178
00180 NumericT operator()(rt_interval<InterfaceType> const & interv,
00181 rt_expr<InterfaceType> const & e,
00182 rt_variable<InterfaceType> const & var) const
00183 {
00184 return quadrature_rule_->eval(interv, e, var);
00185 }
00186
00188 template <id_type id>
00189 NumericT operator()(rt_interval<InterfaceType> const & interv,
00190 rt_expr<InterfaceType> const & e,
00191 ct_variable<id> const & var) const
00192 {
00193 return quadrature_rule_->eval(interv, e, rt_variable<InterfaceType>(id));
00194 }
00195
00196 private:
00197 std::auto_ptr<numerical_quadrature_interface<InterfaceType> > quadrature_rule_;
00198 };
00199 typedef rt_numerical_quadrature<> numerical_quadrature;
00200
00201 }
00202
00203 #endif
00204