test/test_field_integrate.cpp

Go to the documentation of this file.

//
// C++ Implementation: %{MODULE}
//
// Description:
//
//
// Author: %{AUTHOR} <%{EMAIL}>, (C) %{YEAR}
//
// Copyright: See COPYING file that comes with this distribution
//
//

#include <fstream>
#include <tut.h>
#include "myrlog.h"
#include "point.h"
#include "tet.h"
#include "looseoctree.h"
#include "tetmeshbuilder.h"
#include "femaxmesh.h"
#include "LineAtoB.h"
#include "CircleZ.h"
#include "EfieldIntegrand.h"
#include "LoggingFunction.h"
#include "adaptsim.h"
#include "romberg.h"
#include "Epetra_Comm.h"

namespace tut {
    Epetra_Comm& get_comm_world();
}

namespace
{
    using namespace mesh;

    struct test_data
    {
        test_data()
        {
        }
    };
    
    typedef tut::test_group<test_data> testgroup;
    typedef testgroup::object testobject;
    testgroup group("field_integrate");

    mesh::TetMesh* get_mesh4() {
        // Generate Tetmesh
        TetMeshBuilder* builder = new TetMeshBuilder();
        builder->init_coord(6);
        builder->set_coord(0,  1,  0,  0);
        builder->set_coord(1,  0,  1,  0);
        builder->set_coord(2, -1,  0,  0);
        builder->set_coord(3,  0, -1,  0);
        builder->set_coord(4,  0,  0,  1);
        builder->set_coord(5,  0,  0, -1);
        builder->finalize_coord();
        builder->init_tet(4);
        builder->set_tet(0, 0, 1, 3, 4);
        builder->set_tet(1, 1, 2, 3, 4);
        builder->set_tet(2, 0, 1, 3, 5);
        builder->set_tet(3, 1, 2, 3, 5);
        builder->finalize_tet();
        builder->init_bc(8);
        builder->set_bc(0, 1, 4, 0);
        builder->set_bc(0, 1, 5, 0);
        builder->set_bc(0, 3, 4, 0);
        builder->set_bc(0, 3, 5, 0);
        builder->set_bc(1, 2, 4, 0);
        builder->set_bc(1, 2, 5, 0);
        builder->set_bc(2, 3, 4, 0);
        builder->set_bc(2, 3, 5, 0);
        builder->finalize_bc(0);
        TetMesh* mesh = builder->get_mesh();
        delete builder;
        return mesh;
    }

void getdimensions(int nof_cells, mesh::Vector3& box_len, int dim[3]) {
    mesh::Vector3 d;
    d.x = ::pow(static_cast<double>(nof_cells) * box_len.x * box_len.x / box_len.y / box_len.z, 1.0/3.0);
    d.y = d.x * box_len.y / box_len.x;
    d.z = d.x * box_len.z / box_len.x ;

    // make sure that all components are at least 0.5
    for (int i = 0; i < 3; i ++)
        if (d[i] < 0.5)
            d *= 0.5/d[i];

    dim[0] = static_cast<int>(round(d.x));
    dim[1] = static_cast<int>(round(d.y));
    dim[2] = static_cast<int>(round(d.z));
}

        
void export_vtk(const FemaxMesh& femax_mesh, const colarray::Vector<double>& q, const std::string file, int nof_cells) {
    int nofPoints;
    int dim[3];
    mesh::TetMesh* tmesh;
    mesh::Vector3 center;
    mesh::Vector3 box_len;
    mesh::Vector3 value;
    tmesh = &(femax_mesh.get_tet_mesh());

    std::vector<double> xcoords, ycoords, zcoords;
    double curx, cury, curz;
    mesh::Vector3 point;
    mesh::Vector3 ref_point;
    std::vector<mesh::Tet *> tets;
    std::ofstream of;
    int k = 1;

        of.open(file.c_str());
        assert(of.is_open());

        of.precision(6);
        of << "# vtk DataFile Version 2.0" << endl;
        of << "generated using FemaXX" << endl;
        of << "ASCII" << endl << endl;
        of << "DATASET RECTILINEAR_GRID" << endl;

        // get the bounding box from tmesh and calculate the dimensions of the new grid
        tmesh->get_bounding_box(center, box_len);
        ref_point = center - box_len;
        box_len *= 2.0;

        getdimensions(nof_cells, box_len, dim);
        of << "DIMENSIONS " << dim[0] << " " << dim[1] << " " << dim[2] << endl;

        double xstep = box_len.x / dim[0];
        double ystep = box_len.y / dim[1];
        double zstep = box_len.z / dim[2];

        // create three vectors containing the x, y and z coordinates of the points
        // and do the insertion of these into the vtk file
        of << "X_COORDINATES " << dim[0] << " float" << endl;
        for (int i = 0; i <  dim[0]; i++) {
            curx = ref_point.x + (i * xstep);
            xcoords.push_back(curx);
            of << curx << " ";
        }
        of << endl;

        of << "Y_COORDINATES " << dim[1] << " float" << endl;
        for (int i = 0; i <  dim[1]; i++) {
            cury = ref_point.y + (i * ystep);
            ycoords.push_back(cury);
            of << cury << " ";
        }
        of << endl;

        of << "Z_COORDINATES " << dim[2] << " float" << endl;
        for (int i = 0; i <  dim[2]; i++) {
            curz = ref_point.z + (i * zstep);
            zcoords.push_back(curz);
            of << curz << " ";
        }
        of << endl;

        nofPoints = xcoords.size() *  ycoords.size() * zcoords.size() ;
        of << endl << "POINT_DATA " <<  nofPoints <<  endl << endl;

    // writing the electric field values into the vtk file
            of << "VECTORS el_field_mode_" << k << " float" << endl;
            //get the electric field  value for each point
            for (int i = 0; i < dim[2]; i ++) {
                point.z = zcoords.at(i);
                for (int j = 0; j < dim[1]; j ++) {
                    point.y = ycoords.at(j);
                    for (int h = 0; h < dim[0]; h ++) {
                        point.x = xcoords.at(h);
                        value = femax_mesh.get_nedelec_mesh().eval(point, q);
                        of << value.x << " " << value.y << " " << value.z << endl;
                    }
                }
            }
            of << endl;
 
    // writing the magnetic field values into the vtk file
            of << "VECTORS mag_field_mode_" << k << " float" << endl;
            //get the magnetic field  value for each point
            for (int i = 0; i < dim[2]; i ++) {
                point.z = zcoords.at(i);
                for (int j = 0; j < dim[1]; j ++) {
                    point.y = ycoords.at(j);
                    for (int h = 0; h < dim[0]; h ++) {
                        point.x = xcoords.at(h);
                        value = femax_mesh.get_nedelec_mesh().eval_curl(point, q);
                        of << value.x << " " << value.y << " " << value.z << endl;
                    }
                }
            }
            of << endl;
}

    
    template<>
    template<>
    void testobject::test<1>()
    {
        const double tol = 1e-14;
        
        // Generate Tetmesh
        TetMesh* mesh = get_mesh4();
        mesh->log_mesh_info();
        mesh->construct_octree();
        
        // Create FemaxMesh
        FemaxMesh femax_mesh(*mesh, 1, 0x0);
        
        // Fake eigenvector
        colarray::Vector<double> q(femax_mesh.get_nedelec_mesh().get_num_global_mdofs(0));
        ensure(q._n == 1);
        q(0) = 1.0;
        
        // Path with y==0 is orthogonal to E-field
        Vector3 A(0.2, 0.0, -0.5);
        Vector3 B(-0.2, 0.0,  0.5);
        LineAtoB line(A, B);
        EfieldDs<LineAtoB> integrand(line, femax_mesh.get_nedelec_mesh(), q);

        for (int i = 0; i <= 10; ++ i) {
            Vector3 x(line((double)i/10.0));
            Vector3 E(femax_mesh.get_nedelec_mesh().eval(x, q));
            rInfo("E(%f,%f,%f) = (%e,%e,%e)",
                  x.x, x.y, x.z,
                  E.x, E.y, E.z);
            ensure(fabs(E.x) < tol && fabs(E.z) < tol);
        }
      
        int info;
        double gap_voltage = NR::qromb(integrand, 0.0, 1.0, 1e-10, info);
        rInfo("qromb: voltage=%e", gap_voltage);
        ensure(fabs(gap_voltage) < tol);

        gap_voltage = NR::adaptsim(integrand, 0.0, 1.0, 1e-10, info);
        rInfo("adaptsim: voltage=%e, info=%d", gap_voltage, info);
        ensure(info == 0);
        ensure(fabs(gap_voltage) < tol);

        delete mesh;
    }
    
    template<>
    template<>
    void testobject::test<2>()
    {
        const double tol = 1e-14;
        
        // Generate Tetmesh
        TetMesh* mesh = get_mesh4();
        mesh->log_mesh_info();
        mesh->construct_octree();
        
        // Create FemaxMesh
        FemaxMesh femax_mesh(*mesh, 1, 0x0);
        
        // Fake eigenvector
        colarray::Vector<double> q(femax_mesh.get_nedelec_mesh().get_num_global_mdofs(0));
        ensure(q._n == 1);
        q(0) = 1.0;
        
        // Path with y==0 is orthogonal to E-field
        Vector3 A(0.00001, 0.2, -0.5);
        Vector3 B(0.00001, -0.2,  0.5);
        LineAtoB line(A, B);
        EfieldDs<LineAtoB> integrand(line, femax_mesh.get_nedelec_mesh(), q);

        for (int i = 0; i <= 10; ++ i) {
            Vector3 x(line((double)i/10.0));
            Vector3 E(femax_mesh.get_nedelec_mesh().eval(x, q));
            rInfo("E(%f,%f,%f) = (%e,%e,%e)",
                  x.x, x.y, x.z,
                  E.x, E.y, E.z);
        }
      
        int info;
        double gap_voltage_1 = NR::qromb(integrand, 0.0, 1.0, 1e-10, info);
        rInfo("qromb: voltage=%e", gap_voltage_1);

        double gap_voltage_2 = NR::adaptsim(integrand, 0.0, 1.0, 1e-10, info);
        rInfo("adaptsim: voltage=%e, info=%d", gap_voltage_2, info);
        ensure(info == 0);
        ensure(fabs(gap_voltage_2 - gap_voltage_1) < tol);

        delete mesh;
    }

    template<>
    template<>
    void testobject::test<3>()
    {
        const double tol = 1e-14;
        
        // Generate Tetmesh
        TetMesh* mesh = get_mesh4();
        mesh->log_mesh_info();
        mesh->construct_octree();
        
        // Create FemaxMesh
        FemaxMesh femax_mesh(*mesh, 1, 0x0);
        
        // Fake eigenvector
        colarray::Vector<double> q(femax_mesh.get_nedelec_mesh().get_num_global_mdofs(0));
        ensure(q._n == 1);
        q(0) = 100.0;
        
        // Path with y==0 is orthogonal to E-field
        Vector3 A(-0.0005, -0.3, 0.0001);
        Vector3 B(0.0005, 0.5, 0.0001);
        LineAtoB line(A, B);
        EfieldDs<LineAtoB> integrand(line, femax_mesh.get_nedelec_mesh(), q);

        for (int i = 0; i <= 10; ++ i) {
            double t = (double)i/10.0;
            Vector3 x(line(t));
            Vector3 E(femax_mesh.get_nedelec_mesh().eval(x, q));
            rInfo("E(%f,%f,%f), Eds = (%e,%e,%e), %e",
                  x.x, x.y, x.z,
                  E.x, E.y, E.z,
                  integrand(t));
        }
      
        int info;
        double gap_voltage_1 = NR::qromb(integrand, 0.0, 1.0, 1e-15, info);
        rInfo("qromb: voltage=%e", gap_voltage_1);

        double gap_voltage_2 = NR::adaptsim(integrand, 0.0, 1.0, 1e-15, info);
        rInfo("adaptsim: voltage=%e, info=%d", gap_voltage_2, info);
        ensure(info == 0);
        ensure(fabs(gap_voltage_2 - gap_voltage_1) < tol);

        delete mesh;
    }

    template<>
    template<>
    void testobject::test<10>()
    {
        const double int_tol = 2e-8;
        const double a = 0.0;
        //const double b = 4.0*acos(0.0);
        const double b = 6.0;
        
        // Generate Tetmesh
        TetMesh* mesh = get_mesh4();
        mesh->log_mesh_info();
        mesh->construct_octree();
        
        // Create FemaxMesh
        FemaxMesh femax_mesh(*mesh, 1, 0x0);
        
        // Fake eigenvector
        colarray::Vector<double> q(femax_mesh.get_nedelec_mesh().get_num_global_mdofs(0));
        ensure(q._n == 1);
        q(0) = 1.0;
        
        // Circular path
        CircleZ curve(0.5, 0.01);
        EfieldDs<CircleZ> integrand(curve, femax_mesh.get_nedelec_mesh(), q);
        LoggingFunction<EfieldDs<CircleZ>, double> my_integrand(integrand);

        const int max_steps = 100;
        for (int i = 0; i <= max_steps; ++ i) {
            double t = (b-a) * (double)i/max_steps + a;
            Vector3 x(curve(t));
            Vector3 E(femax_mesh.get_nedelec_mesh().eval(x, q));
            rInfo("E(%f,%f,%f), Eds = (%e,%e,%e), %e",
                  x.x, x.y, x.z,
                  E.x, E.y, E.z,
                  integrand(t));
        }
      
        int info;
        double gap_voltage_1 = NR::qromb(my_integrand, a, b, int_tol, info);
        rInfo("qromb: voltage=%e", gap_voltage_1);
        rInfo("num_eval=%d, y_min=%e, y_max=%e", 
              my_integrand.get_num_eval(), 
              my_integrand.get_y_min(), 
              my_integrand.get_y_max());
        ensure(info == 0);

        my_integrand.reset();
        if (tut::get_comm_world().MyPID() == 0)
            my_integrand.enable_plot("adaptsim.tab");
        double gap_voltage_2 = NR::adaptsim(my_integrand, a, b, int_tol, info);
        rInfo("adaptsim: voltage=%e, info=%d", gap_voltage_2, info);
        rInfo("num_eval=%d, y_min=%e, y_max=%e", 
              my_integrand.get_num_eval(), 
              my_integrand.get_y_min(), 
              my_integrand.get_y_max());
        ensure(info == 0);
        
        my_integrand.reset();
        double gap_voltage_3 = NR::qtrap(my_integrand, a, b, int_tol, info);
        rInfo("qtrap: voltage=%e, info=%d", gap_voltage_3, info);
        rInfo("num_eval=%d, y_min=%e, y_max=%e", 
              my_integrand.get_num_eval(), 
              my_integrand.get_y_min(), 
              my_integrand.get_y_max());
        ensure(info == 0);
        
        if (!(fabs(gap_voltage_2 - gap_voltage_1) < 1e3*int_tol))
            rWarning("Results from qromb and adaptsim differ. qromb returns too early.");
        // ensure(fabs(gap_voltage_2 - gap_voltage_1) < 1e3*int_tol);
        ensure(fabs(gap_voltage_3 - gap_voltage_2) < 1e3*int_tol);

        delete mesh;
    }    
    
    template<>
    template<>
    void testobject::test<-1>()
    {
        // Generate Tetmesh
        TetMesh* mesh = get_mesh4();
        mesh->log_mesh_info();
        mesh->construct_octree();
        
        // Create FemaxMesh
        FemaxMesh femax_mesh(*mesh, 1, 0x0);
        
        // Fake eigenvector
        colarray::Vector<double> q(femax_mesh.get_nedelec_mesh().get_num_global_mdofs(0));
        ensure(q._n == 1);
        q(0) = 1.0;
        
        // Export mesh/field to vtk file
        if (tut::get_comm_world().MyPID() == 0)
            export_vtk(femax_mesh, q, "test_fieldintegrate.vtk", 1000);
        
        delete mesh;
    }
    
    template<>
    template<>
    void testobject::test<-2>()
    {
        // Generate Tetmesh
        TetMesh* mesh = get_mesh4();
        mesh->log_mesh_info();
        mesh->construct_octree();
        
        // Create FemaxMesh
        FemaxMesh femax_mesh(*mesh, 1, 0x0);
        
        // Fake eigenvector
        colarray::Vector<double> q(femax_mesh.get_nedelec_mesh().get_num_global_mdofs(0));
        ensure(q._n == 1);
        q(0) = 1.0;
        
        // Circular path
        CircleZ curve(0.5, 0.01);
        EfieldDs<CircleZ> integrand(curve, femax_mesh.get_nedelec_mesh(), q);

        for (double int_tol = 1.0; int_tol >= 1e-15; int_tol /= 10.0) {
            int info;
            double gap_voltage = NR::qromb(integrand, 0.0, 6.0, int_tol, info);
            rWarning("qromb: tol=%e, voltage=%.12e, info=%d", int_tol, gap_voltage, info);
            ensure(info == 0);
        }
        
        for (double int_tol = 1.0; int_tol >= 1e-15; int_tol /= 10.0) {
            int info;
            double gap_voltage = NR::qtrap(integrand, 0.0, 6.0, int_tol, info);
            rWarning("qtrap: tol=%e, voltage=%.12e, info=%d", int_tol, gap_voltage, info);
            ensure(info == 0);
        }
        
        for (double int_tol = 1.0; int_tol >= 1e-15; int_tol /= 10.0) {
            int info;
            double gap_voltage = NR::adaptsim(integrand, 0.0, 6.0, int_tol, info);
            rWarning("adaptsim: tol=%e, voltage=%.12e, info=%d", int_tol, gap_voltage, info);
            ensure(info == 0);
        }
        delete mesh;
    }

}

---