test/test_line_integral.cpp

Go to the documentation of this file.
00001 //
00002 // C++ Implementation: test_romberg
00003 //
00004 // Description: unit tests for line integrals
00005 //
00006 //
00007 // Author: Roman Geus <geus@maxwell>, (C) 2004
00008 //
00009 // Copyright: See COPYING file that comes with this distribution
00010 //
00011 //
00012 
00013 #include <cmath>
00014 #include <tut.h>
00015 #include "CircleZ.h"
00016 #include "LineAtoB.h"
00017 #include "romberg.h"
00018 
00019 using mesh::Vector3;
00020 
00021 template <typename Curve>
00022 struct Curve_length {
00023     explicit Curve_length(Curve curve)
00024         : curve_(curve)
00025     {
00026     }
00027     double operator()(double t) {
00028         Vector3 d(curve_.deriv(t));
00029         return sqrt(d.x*d.x + d.y*d.y + d.z*d.z);
00030     }
00031     Curve curve_;
00032 };
00033 
00036 template <typename Curve>
00037 class MyEfieldDs {
00038 public:
00046     MyEfieldDs(Curve curve, Vector3 E)
00047         : curve_(curve),
00048           E_(E)
00049     { }
00051     double operator()(const double t) {
00052         Vector3 field(E_);
00053         return field.dot_product(curve_.deriv(t));
00054     }
00055 private:
00057     Curve curve_;
00059     const Vector3 E_;
00060 };
00061 
00062 namespace
00063 {
00064 
00065     struct testdata
00066     {
00067         testdata()
00068             : pi(std::acos(-1.0)),
00069               tol(1e-10)
00070         {
00071         }
00072         const double pi;
00073         const double tol;
00074     };
00075 
00076     typedef tut::test_group<testdata> testgroup;
00077     typedef testgroup::object testobject;
00078     testgroup romberg_group("line_integral");
00079 
00081     template<>
00082     template<>
00083     void testobject::test<1>()
00084     {
00085         int info;
00086         const double radius = 2.0;
00087         CircleZ curve(radius, 0.0);
00088         Curve_length<CircleZ> f_len(curve);
00089         double I = NR::qromb<Curve_length<CircleZ> >(f_len, 0.0, 2.0*pi, 1e-12, info);
00090         double expected = 2.0*pi*radius;
00091         ensure(fabs(I - expected) <= tol);
00092     }
00093 
00095     template<>
00096     template<>
00097     void testobject::test<2>()
00098     {
00099         int info;
00100         Vector3 a(2.3, 7.8, -1.2);
00101         Vector3 b(-2.0, 5.1, 3.7);
00102         Vector3 e(-2.5, -8.9, 4.5);
00103         LineAtoB curve(a, b);
00104         MyEfieldDs<LineAtoB> e_field_integrand(curve, e);
00105         double I = NR::qromb<MyEfieldDs<LineAtoB> >(e_field_integrand, 0.0, 1.0, 1e-12, info);
00106         Vector3 d(b-a);
00107         double Eproj = e.dot_product(d) / d.length();
00108         double expected = Eproj*d.length();
00109         ensure(fabs(I - expected) <= tol);
00110     }
00111 
00112 };

Generated on Fri Oct 26 13:35:13 2007 for FEMAXX (Finite Element Maxwell Eigensolver) by  doxygen 1.4.7