test/test_line_integral.cpp

Go to the documentation of this file.

//
// C++ Implementation: test_romberg
//
// Description: unit tests for line integrals
//
//
// Author: Roman Geus <geus@maxwell>, (C) 2004
//
// Copyright: See COPYING file that comes with this distribution
//
//

#include <cmath>
#include <tut.h>
#include "CircleZ.h"
#include "LineAtoB.h"
#include "romberg.h"

using mesh::Vector3;

template <typename Curve>
struct Curve_length {
    explicit Curve_length(Curve curve)
        : curve_(curve)
    {
    }
    double operator()(double t) {
        Vector3 d(curve_.deriv(t));
        return sqrt(d.x*d.x + d.y*d.y + d.z*d.z);
    }
    Curve curve_;
};

template <typename Curve>
class MyEfieldDs {
public:
    MyEfieldDs(Curve curve, Vector3 E)
        : curve_(curve),
          E_(E)
    { }
    double operator()(const double t) {
        Vector3 field(E_);
        return field.dot_product(curve_.deriv(t));
    }
private:
    Curve curve_;
    const Vector3 E_;
};

namespace
{

    struct testdata
    {
        testdata()
            : pi(std::acos(-1.0)),
              tol(1e-10)
        {
        }
        const double pi;
        const double tol;
    };

    typedef tut::test_group<testdata> testgroup;
    typedef testgroup::object testobject;
    testgroup romberg_group("line_integral");

    template<>
    template<>
    void testobject::test<1>()
    {
        int info;
        const double radius = 2.0;
        CircleZ curve(radius, 0.0);
        Curve_length<CircleZ> f_len(curve);
        double I = NR::qromb<Curve_length<CircleZ> >(f_len, 0.0, 2.0*pi, 1e-12, info);
        double expected = 2.0*pi*radius;
        ensure(fabs(I - expected) <= tol);
    }

    template<>
    template<>
    void testobject::test<2>()
    {
        int info;
        Vector3 a(2.3, 7.8, -1.2);
        Vector3 b(-2.0, 5.1, 3.7);
        Vector3 e(-2.5, -8.9, 4.5);
        LineAtoB curve(a, b);
        MyEfieldDs<LineAtoB> e_field_integrand(curve, e);
        double I = NR::qromb<MyEfieldDs<LineAtoB> >(e_field_integrand, 0.0, 1.0, 1e-12, info);
        Vector3 d(b-a);
        double Eproj = e.dot_product(d) / d.length();
        double expected = Eproj*d.length();
        ensure(fabs(I - expected) <= tol);
    }

};

---