test/test_field_interp.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 <stdlib.h>
#include <iostream>
#include <fstream>
#include <sstream>
#include <unistd.h>
#include <tut.h>
#include "myrlog.h"
#include "tetmesh/point.h"
#include "tetmesh/tet.h"
#include "colarray/vector.h"
#include "looseoctree.h"
#include "tetmeshbuilder.h"
#include "broadcastmeshbuilder.h"
#include "netgenreader.h"
#include "femaxmesh.h"
#include "LinearFieldInterpolation.h"

namespace tut {
    Epetra_Comm& get_comm_world();
}

namespace
{
    using namespace mesh;
    using namespace postprocess;

    struct test_data
    {
        test_data()
        {
        }
        static mesh::TetMesh* tmesh_;
        static FemaxMesh* femax_mesh_;
    };

    mesh::TetMesh* test_data::tmesh_ = 0;
    FemaxMesh* test_data::femax_mesh_ = 0;
    
    typedef tut::test_group<test_data> testgroup;
    typedef testgroup::object testobject;
    testgroup group("field_interp");

    template<>
    template<>
    void testobject::test<1>()
    {
        // Read mesh data structure
        std::string mesh_file_name = "test/octree/cop40k.ng";
        rAssert(!access(mesh_file_name.c_str(), R_OK));

        std::ifstream istr(mesh_file_name.c_str());

        Epetra_Comm& comm(tut::get_comm_world());
        TetMeshBuilder* builder = new TetMeshBuilder();
#ifdef HAVE_MPI
        BroadcastMeshBuilder* broadcast = new BroadcastMeshBuilder(&comm, builder);
        NetgenReader* reader(0);
        if (comm.MyPID() == 0) {
            reader = new NetgenReader(&istr, broadcast);
            reader->read();
        } else
            broadcast->receiver();
        delete broadcast;
        if (comm.MyPID() == 0)
            delete reader;
#else
        NetgenReader* reader = new NetgenReader(&istr, builder);
        reader->read();
        delete reader;
#endif
        tmesh_ = builder->get_mesh();
        delete builder;
        tmesh_->log_mesh_info();

        // Construct Octree
        tmesh_->construct_octree();
        
        // Construct FemaxMesh
        const int element_order = 1;
        const int sym_plane_config = 6;
        femax_mesh_ = new FemaxMesh(*tmesh_, element_order, sym_plane_config);
    }

    template<>
    template<>
    void testobject::test<10>()
    {
        const double diff_tol = 1e-14;
    
        // Create fake random eigenvector
        colarray::Vector<double> q(femax_mesh_->get_nedelec_mesh().get_num_global_mdofs(0));
        for (int i = 0; i < q._n; ++ i)
            q(i) = static_cast<double>(random()) / RAND_MAX;
        
        // Create LinearFieldInterpolation
        LinearFieldInterpolation eval_interp(femax_mesh_->get_nedelec_mesh(), q);
        
        // Evaluate and compare field at mesh points    
        for (TetMesh::point_iterator it = tmesh_->point_begin(); it != tmesh_->point_end(); ++ it) {
            Point& point = *it;
            Vector3 f1 = eval_interp.eval(point.get_coord());
            Vector3 f2 = femax_mesh_->get_nedelec_mesh().eval(point.get_coord(), q);
            Vector3 d = f1 - f2;
            double diff = d.length();
            ensure(diff < diff_tol);
        }
    }

    template<>
    template<>
    void testobject::test<11>()
    {
        const double diff_tol = 1e-12;
    
        // Create fake random eigenvector
        colarray::Vector<double> q(femax_mesh_->get_nedelec_mesh().get_num_global_mdofs(0));
        for (int i = 0; i < q._n; ++ i)
            q(i) = static_cast<double>(random()) / RAND_MAX;
        
        // Create LinearFieldInterpolation
        LinearFieldInterpolation eval_interp(femax_mesh_->get_nedelec_mesh(), q);
        
        // Evaluate and compare field at mesh points    
        for (TetMesh::tet_iterator it = tmesh_->tet_begin(); it != tmesh_->tet_end(); ++ it) {
            Tet* tet = &*it;
            // Interpolate field at mid-point of Tet
            Vector3 mid_coord(Vector3::ZERO);
            for (int i = 0; i < 4; ++ i)
                mid_coord += tet->get_corner(i)->get_coord();
            mid_coord /= 4.0;
            Vector3 f1 = eval_interp.eval(mid_coord);
            // Build average of field evaluated at the four corners
            Vector3 f2(Vector3::ZERO);
            for (int i = 0; i < 4; ++ i)
                f2 += eval_interp.eval(tet->get_corner(i)->get_coord());
            f2 /= 4.0;
            // Compare
            Vector3 d = f1 - f2;
            double diff = d.length();
            ensure(diff < diff_tol);
        }
    }

    template<>
    template<>
    void testobject::test<50>()
    {
        delete femax_mesh_;
        delete tmesh_;
    }

}

---