test/test_octree.cpp

Go to the documentation of this file.

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

#include <iostream>
#include <fstream>
#include <sstream>
#include <unistd.h>
#include <tut.h>
#include "myrlog.h"
#include "point.h"
#include "tet.h"
#include "looseoctree.h"
#include "tetmeshbuilder.h"
#include "broadcastmeshbuilder.h"
#include "netgenreader.h"

namespace tut {
    Epetra_Comm& get_comm_world();
}

namespace
{
    using namespace mesh;

    struct test_data
    {
        test_data()
        {
        }
        static mesh::TetMesh* tmesh;
    };

    mesh::TetMesh* test_data::tmesh = 0;

    typedef tut::test_group<test_data> testgroup;
    typedef testgroup::object testobject;
    testgroup group("octree");

    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();
    }

    template<>
    template<>
    void testobject::test<10>()
    {
        // test if all mesh points are inside mesh
        for (mesh::TetMesh::point_iterator it = tmesh->point_begin(); it != tmesh->point_end(); ++ it) {
            mesh::Vector3 p = (*it).get_coord();
            bool inside = tmesh->is_inside(p);
            if (!inside)
                rError("point %d outside.", (*it).get_id());
            ensure(inside);
        }
    }

    template<>
    template<>
    void testobject::test<11>()
    {
        LooseOctree<Tet>* octree = tmesh->get_octree();

        // Compute number of tets connected to each node
        vector<int> conn_tet(tmesh->get_nof_points());
        std::fill(conn_tet.begin(), conn_tet.end(), 0);
        for (TetMesh::tet_iterator it = tmesh->tet_begin(); it != tmesh->tet_end(); ++ it) {
            for (int i=0; i < 4; ++ i)
                conn_tet[(*it).get_corner_id(i)] ++;
        }

        for (TetMesh::point_iterator it = tmesh->point_begin(); it != tmesh->point_end(); ++ it) {
            Vector3 p = (*it).get_coord();

            // Compute list of tets potentially containing p
            std::vector<Tet *> nodes;
            int nof_visited = 0;
            octree->find_nodes_by_point(p, nodes, nof_visited, &Tet::is_inside_bounding_box_gracious);


            // Compute list of tets containing point p
            std::vector<Tet*> tets;
            for (std::vector<Tet *>::iterator it2 = nodes.begin(); it2 != nodes.end(); ++ it2) {
                Tet *tet = dynamic_cast<Tet*>(*it2);
                if (tet->is_inside_gracious(p))
                    tets.push_back(tet);
            }

            int nof_found = static_cast<int>(tets.size());
            if (nof_found != conn_tet[(*it).get_id()])
                rError("point %d: found %d tets, expected %d.",
                       (*it).get_id(),
                       nof_found,
                       conn_tet[(*it).get_id()]);
            ensure(nof_found == conn_tet[(*it).get_id()]);
        }
    }

    template<>
    template<>
    void testobject::test<20>()
    {
        for (TetMesh::tet_iterator it = tmesh->tet_begin(); it != tmesh->tet_end(); ++ it) {
            Tet* tet(&*it);
            for (int i=0; i < 4; ++ i)
                ensure(tet->is_inside(tet->get_corner(i)->get_coord()));
        }
    }

    template<>
    template<>
    void testobject::test<21>()
    {
        for (TetMesh::tet_iterator it = tmesh->tet_begin(); it != tmesh->tet_end(); ++ it) {
            Tet* tet(&*it);
            for (int i=0; i < 4; ++ i)
                ensure(tet->is_inside_gracious(tet->get_corner(i)->get_coord()));
        }
    }

    template<>
    template<>
    void testobject::test<22>()
    {
        for (TetMesh::tet_iterator it = tmesh->tet_begin(); it != tmesh->tet_end(); ++ it) {
            Tet* tet(&*it);
            for (int i=0; i < 4; ++ i)
                ensure(tet->is_inside_bounding_box(tet->get_corner(i)->get_coord()));
        }
    }

    template<>
    template<>
    void testobject::test<23>()
    {
        for (TetMesh::tet_iterator it = tmesh->tet_begin(); it != tmesh->tet_end(); ++ it) {
            Tet* tet(&*it);
            for (int i=0; i < 4; ++ i)
                ensure(tet->is_inside_bounding_box_gracious(tet->get_corner(i)->get_coord()));
        }
    }

    template<>
    template<>
    void testobject::test<50>()
    {
        // delete mesh data structure (including octree)
        delete tmesh;
    }

}

---