src/tetmesh/tetmesh.cpp

Go to the documentation of this file.

/***************************************************************************
                          tetmesh.cpp  -  description
                             -------------------
    begin                : Mon Dec 15 2003
    copyright            : (C) 2003 by Roman Geus
    email                : roman.geus@psi.ch
***************************************************************************/

/***************************************************************************
 *                                                                         *
 *   This program is free software; you can redistribute it and/or modify  *
 *   it under the terms of the GNU General Public License as published by  *
 *   the Free Software Foundation; either version 2 of the License, or     *
 *   (at your option) any later version.                                   *
 *                                                                         *
 ***************************************************************************/

#include <cassert>
#include <iomanip>
#include <new>
#include <vector>
#include <stdexcept>
#include "myrlog.h"
#include "pbedefs.h"
#include "tetmesh.h"

namespace mesh {

using namespace std;

TetMesh* TetMesh::instance_ = 0;        // Pointer to TetMesh singleton.

TetMesh::TetMesh()
    : _nofTet(0),
      _nofPoint(0),
      _nofEdge(0),
      _nofFace(0),
      _nofSym(0),
      _tetCounter(0),
      _octree(0)
{
    assert(instance_ == 0);
    instance_ = this;
}

TetMesh::~TetMesh() {
    // delete octree
    if (_octree)
        delete _octree;
    instance_ = 0;
}


void TetMesh::load_materials(const char* fileName){
  _materials = new Materials(fileName);
}

Materials* TetMesh::get_materials() const{
  return _materials;
}


void TetMesh::initPoint(int nofPoint) {
    _nofPoint = nofPoint;
    _points.reserve(nofPoint);
}

void TetMesh::initTet(int nofTet) {
    _nofTet = nofTet;
    _tets.reserve(nofTet);
}

void TetMesh::insertPoint(int id, double x, double y, double z) {
    assert(static_cast<size_t>(id) == _points.size());
    _points.push_back(Point(id, x, y, z));
}

  void TetMesh::insertTet(id_t global_id, id_t id0, id_t id1, id_t id2, id_t id3, id_t material) {

    id_t id = _tetCounter ++;
    id_t corners[4] = {id0, id1, id2, id3};
    assert(id == _tets.size());
    _tets.push_back(Tet(*this, id, corners, material));
}

void TetMesh::generateEdges() {
    const int edgeNodeLocal[6][2] = {{0, 1}, {0, 2}, {0, 3},
                                     {1, 2}, {1, 3}, {2, 3}};
    id_t node0_id;
    id_t node1_id;
    id_t edge_id;

    for (tet_iterator iter = tet_begin(); iter != tet_end(); ++ iter) {
        for (int j = 0; j < 6; j ++) {
            node0_id = (*iter).get_corner_id(edgeNodeLocal[j][0]);
            node1_id = (*iter).get_corner_id(edgeNodeLocal[j][1]);
            edge_id = insertEdge(node0_id, node1_id);
            (*iter).set_edge_id(j, edge_id);
            // FIXME: orientation edges must be consistent globally
            (*iter).set_edge_orientation(j, get_edge(edge_id)->get_start_id() == node0_id);
        }
    }
}

void TetMesh::generateFaces() {
    int node0, node1, node2;
    for (tet_iterator iter = tet_begin(); iter != tet_end(); ++ iter) {
        for (int j = 0; j < 4; j ++) {
            Tet::get_face_point_lids(j, node0, node1, node2);
            id_t face_id = insertFace((*iter).get_id(),
                                      (*iter).get_corner_id(node0),
                                      (*iter).get_corner_id(node1),
                                      (*iter).get_corner_id(node2));
            (*iter).set_face_id(j, face_id);
        }
    }
    // dump_face_map(cerr);
}

id_t TetMesh::insertEdge(id_t node0, id_t node1) {
    if (node0 > node1)    /* force n0 < n1 */
        std::swap<id_t>(node0, node1);
    EdgeMapKey key(node0, node1);
    EdgeMap::iterator iter = _edgeMap.find(key);
    if (iter == _edgeMap.end()) {
        _edges.push_back(Edge(_nofEdge, node0, node1));
        id_t edge = _nofEdge;
        _edgeMap[key] = edge;
        _nofEdge += 1;
        return edge;
    } else
        return (*iter).second;
}

id_t TetMesh::insertFace(id_t tet, id_t node0, id_t node1, id_t node2) {
    // sort node ids in ascending order
    if (node0 > node1)
        std::swap<id_t>(node0, node1);
    if (node1 > node2)
        std::swap<id_t>(node1, node2);
    if (node0 > node1)
        std::swap<id_t>(node0, node1);

    // generate key for map
    FaceMapKey key(node0, node1, node2);
    // find corresponding entry in map
    FaceMap::iterator iter = _faceMap.find(key);
    
    //cout << "TetMesh inserting "<<tet<<" local nodes "<<node0<<","<<node1<<","<<node2<<endl;
    
    if (iter == _faceMap.end()) {
        // no entry found: create new face and insert it in map
        _faces.push_back(Face(_nofFace, node0, node1, node2, tet));
        _faceMap[key] = _nofFace;
        return _nofFace ++;
    } else {
        // set 2nd parent tet for face
        id_t face_id = (*iter).second;
        get_face(face_id)->set_tet_id(tet);
        // return id of face stored in map
        return face_id;
    }
}

id_t TetMesh::lookupEdge(id_t node0, id_t node1) {
    if (node0 > node1)    /* force n0 < n1 */
        std::swap<id_t>(node0, node1);

    EdgeMap::iterator iter = _edgeMap.find(EdgeMapKey(node0, node1));
    if (iter == _edgeMap.end())
        throw std::runtime_error("Accessed non-existing edge.");
    return (*iter).second;
}

id_t TetMesh::lookupFace(id_t node0, id_t node1, id_t node2) {
    // sort node ids in ascending order
    if (node0 > node1)
        std::swap<id_t>(node0, node1);
    if (node1 > node2)
        std::swap<id_t>(node1, node2);
    if (node0 > node1)
        std::swap<id_t>(node0, node1);

    FaceMap::iterator iter = _faceMap.find(FaceMapKey(node0, node1, node2));
    if (iter == _faceMap.end())    
        throw std::runtime_error("Accessed non-existing face.");
    return (*iter).second;
}

void TetMesh::dump_face_map(ostream& os) {
    os << "face map dump:" << endl;
    FaceMap::iterator it(_faceMap.begin());
    for (; it != _faceMap.end(); ++ it) {
        os << (*it).first.first << ", " <<  
              (*it).first.second << ", " <<  
              (*it).first.third << ": " <<  
              (*it).second <<  
              endl;
    }
}

void TetMesh::setBoundaryOnFace(id_t face_id) {
    Face* face = get_face(face_id);
    id_t corners[3] = { face->get_corner_id(0),
                        face->get_corner_id(1),
                        face->get_corner_id(2) };
    face->set_boundary(true);

    for (int i = 0; i < 3; i ++)
        get_point(corners[i])->set_boundary(true);

    get_edge(lookupEdge(corners[0], corners[1]))->set_boundary(true);
    get_edge(lookupEdge(corners[0], corners[2]))->set_boundary(true);
    get_edge(lookupEdge(corners[1], corners[2]))->set_boundary(true);
}

void TetMesh::setSymmetryPlaneOnFace(id_t face_id, int sym_id) {
    Face* face = get_face(face_id);
    id_t corners[3] = { face->get_corner_id(0),
                        face->get_corner_id(1),
                        face->get_corner_id(2) };
    face->set_symmetry_plane(sym_id);

    for (int i = 0; i < 3; i ++)
        get_point(corners[i])->set_symmetry_plane(sym_id);

    get_edge(lookupEdge(corners[0], corners[1]))->set_symmetry_plane(sym_id);
    get_edge(lookupEdge(corners[0], corners[2]))->set_symmetry_plane(sym_id);
    get_edge(lookupEdge(corners[1], corners[2]))->set_symmetry_plane(sym_id);
}

void TetMesh::finalize_mesh(int nofSym) {
    _nofSym = nofSym;
    // Edge and face maps are no longer needed. Erase all their data.
    // FIXME: Hopefully clear() really frees up the allocated memory (at least glibc++ 
    //        delallocates all memory).
    _edgeMap.clear();
    _faceMap.clear();
}

Tet* TetMesh::get_tet(id_t id) const { 
    return const_cast<Tet *>(&_tets[id]); 
}

Face* TetMesh::get_face(id_t id) const {
    return const_cast<Face *>(&_faces[id]); 
}

Edge* TetMesh::get_edge(id_t id) const { 
    return const_cast<Edge *>(&_edges[id]); 
}

Point* TetMesh::get_point(id_t id) const { 
    return const_cast<Point *>(&_points[id]);
}

void TetMesh::log_mesh_info() {
    // estimate mem consumption of tet, face, edge and point objects
    unsigned long size = _nofTet * sizeof(Tet) * 11 / 10;
    size += _nofFace * sizeof(Face) * 11 / 10;
    size += _nofEdge * sizeof(Edge) * 11 / 10;
    size += _nofPoint * sizeof(Point) * 11 / 10;
    // estimate mem consumption of edge and face maps (raw data + 10% overhead)
    size += _nofEdge * sizeof(EdgeMapKey) + sizeof(id_t) * 11 / 10;
    size += _nofFace * sizeof(FaceMapKey) + sizeof(id_t) * 11 / 10;
    rInfo("Mesh stats: #tets=%d, #faces=%d, #edges=%d, #nodes=%d, #symplanes=%d, mem=%s",
          _nofTet, _nofFace, _nofEdge, _nofPoint, _nofSym, rlog::bytes2str(size).c_str());
    rDebug("sizeof(Tet)=%d, sizeof(Face)=%d, sizeof(Edge)=%d, sizeof(Point)=%d",
           static_cast<int>(sizeof(Tet)),
           static_cast<int>(sizeof(Face)),
           static_cast<int>(sizeof(Edge)),
           static_cast<int>(sizeof(Point)));
}

void TetMesh::get_bounding_box(Vector3& center, Vector3& extent) const {
    Vector3 min, max, coord;

    Point* point = get_point(0);
    min = point->get_coord();
    max = min;
    for (id_t i = 1; i < get_nof_points(); i ++) {
        coord = get_point(i)->get_coord();
        min.make_floor(coord);
        max.make_ceil(coord);
    }
    center = max.mid_point(min);
    extent = 0.5 * (max - min);
}

void TetMesh::construct_octree() {
    rAssert(_octree == 0);
    
    // bounding box of whole mesh
    Vector3 center, extent;
    get_bounding_box(center, extent);

    // octree length is the largest component of extent
    double length = extent.x;
    if (extent.y > length)
        length = extent.y;
    if (extent.z > length)
        length = extent.z;
  
    // construct LooseOctree with large max depth
    _octree = new LooseOctree<Tet>(center, length, 100);

    // loop over all tets and insert them in the octree
    for (id_t t = 0; t < get_nof_tets(); t ++) {
        _octree->add_node(get_tet(t));
    }
    int max_depth, nof_octants, nof_nodes;
    double avg_node_depth;
    _octree->get_statistics(max_depth, avg_node_depth, nof_octants, nof_nodes);
    rInfo("Octree stats: max. depth=%d, #octants=%d, #nodes=%d, avg. #node per octant=%.2f, avg. node depth=%.2f",
          max_depth, nof_octants, nof_nodes, static_cast<double>(nof_nodes)/nof_octants, avg_node_depth);

#if 0
    // Vector3 p (0.000001, 0.000001, 0.000001);
    // Vector3 p (0.000000, 0.000000, 0.000000);
    // Vector3 p (1.00000071059, 0.500001458653, 0.106184447056);
    Vector3 p (0.1, 0.0, 2.6);
    std::vector<OctreeNode *> nodes;
    int nof_visited;

    cout << "Finding tets containing " << p << "\n\n";

    tet_iterator tet_it;
    cout << "Exhausive search yields (is_inside_bounding_box_gracious)\n";
    for (tet_it = tet_begin(); tet_it != tet_end(); ++ tet_it) {
        Tet* tet = *tet_it;
        if (tet->is_inside_bounding_box_gracious(p)) {
            std::cout << "inside tet " << tet->get_id() << "\n";
            for (int i = 0; i < 4; i ++) {
                cout << "   " << tet->get_corner(i)->get_coord() << "\n";
            }
        }
    }

    cout << "Exhausive search yields (is_inside_bounding_box)\n";
    for (tet_it = tet_begin(); tet_it != tet_end(); ++ tet_it) {
        Tet* tet = *tet_it;
        if (tet->is_inside_bounding_box(p)) {
            std::cout << "inside tet " << tet->get_id() << "\n";
            for (int i = 0; i < 4; i ++) {
                cout << "   " << tet->get_corner(i)->get_coord() << "\n";
            }
        }
    }

    cout << "Exhausive search yields (is_inside)\n";
    for (tet_it = tet_begin(); tet_it != tet_end(); ++ tet_it) {
        Tet* tet = *tet_it;
        if (tet->is_inside(p)) {
            std::cout << "inside tet " << tet->get_id() << "\n";
            for (int i = 0; i < 4; i ++) {
                cout << "   " << tet->get_corner(i)->get_coord() << "\n";
            }
        }
    }

    cout << "\nOctree: is_inside_bounding_box\n";
    nof_visited = 0; nodes.resize(0);
    _octree->find_nodes_by_point(p, nodes, nof_visited,
                                 (&OctreeNode::is_inside_bounding_box));
    std::cout << "nof. nodes visited=" << nof_visited << "\n";
    std::cout << "nof. nodes found=" << nodes.size() << "\n";
    std::vector<OctreeNode *>::iterator it;
    for (it = nodes.begin(); it != nodes.end(); it ++) {
        Tet *tet = dynamic_cast<Tet*>(*it);
        if (tet->is_inside(p)) {
            std::cout << "inside tet " << tet->get_id() << "\n";
            for (int i = 0; i < 4; i ++) {
                cout << "   " << tet->get_corner(i)->get_coord() << "\n";
            }
        }
    }

    cout << "\nOctree: is_inside_bounding_box_gracious\n";
    nof_visited = 0; nodes.resize(0);
    _octree->find_nodes_by_point(p, nodes, nof_visited,
                                 (&OctreeNode::is_inside_bounding_box_gracious));
    std::cout << "nof. nodes visited=" << nof_visited << "\n";
    std::cout << "nof. nodes found=" << nodes.size() << "\n";
    for (it = nodes.begin(); it != nodes.end(); it ++) {
        Tet *tet = dynamic_cast<Tet*>(*it);
        if (tet->is_inside(p)) {
            std::cout << "inside tet " << tet->get_id() << "\n";
            for (int i = 0; i < 4; i ++) {
                cout << "   " << tet->get_corner(i)->get_coord() << "\n";
            }
        }
    }
   
    cout << "\nOctree: Tet::is_inside\n";
    nof_visited = 0; nodes.resize(0);
    _octree->find_nodes_by_point(p, nodes, nof_visited,
                                 (LooseOctree::filter)(&Tet::is_inside));
    std::cout << "nof. nodes visited=" << nof_visited << "\n";
    std::cout << "nof. nodes found=" << nodes.size() << "\n";
    for (it = nodes.begin(); it != nodes.end(); it ++) {
        Tet *tet = dynamic_cast<Tet*>(*it);
        if (tet->is_inside(p)) {
            std::cout << "inside tet " << tet->get_id() << "\n";
            for (int i = 0; i < 4; i ++) {
                cout << "   " << tet->get_corner(i)->get_coord() << "\n";
            }
        }
    }

    cout << "\nOctree, find_one: Tet::is_inside\n";
    nof_visited = 0;
    // this is a dangerous typecast: it's only valid since we know all nodes in the
    // octree are Tets.
    OctreeNode* node = _octree->find_one_node_by_point(p, nof_visited, (LooseOctree::filter)(&Tet::is_inside));
    std::cout << "nof. nodes visited=" << nof_visited << "\n";
    if (node) {
        std::cout << "nof. nodes found=" << 1 << "\n";
        Vector4 p_simplex;
        dynamic_cast<Tet *>(node)->cartesian_to_simplex(p, p_simplex);
        cout << "is_inside(p) is " << dynamic_cast<Tet *>(node)->is_inside(p) << "\n";
        std::cout << "inside tet " << dynamic_cast<Tet *>(node)->get_id() << "\n";
        for (int i = 0; i < 4; i ++) {
            cout << "   " << dynamic_cast<Tet*>(node)->get_corner(i)->get_coord() << "\n";
        }
        cout << "p_simplex=" << p_simplex << "\n";
    }

    {
        Vector3 center, extent;
        get_tet(6488)->get_bounding_box(center, extent);
        cout.precision(17);
        cout << "\n\ncenter: " << center << "\nextent: " << extent << "\n";
        center -= p;
        center.x = fabs(center.x);
        center.y = fabs(center.y);
        center.z = fabs(center.z);
        cout << "|c-p|: " << center
             << "\nextend - |c-p|:" << extent - center
             << "\nextend_ - |c-p|:" << extent*(1+1e-13) - center
             << "\nbool: " << (center <= extent) << "\n";
        _octree->add_node(get_tet(6488))->dump(cout);
        cout << "eps_mach=" << std::numeric_limits<double>::epsilon() << "\n";
    }
#endif
}

void TetMesh::find_tets_by_point(const Vector3& p, vector<Tet *>& tets) const {
    assert(_octree);

    pbe_start(132, "Octree traversal");

    // Get list of nodes (Tets) whose bounding box contain point p
    std::vector<Tet*> nodes;
    int nof_visited = 0;
    _octree->find_nodes_by_point(p, nodes, nof_visited, &Tet::is_inside_bounding_box_gracious);

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

    pbe_stop(132);
    // rDebug("find_tets_by_point: nof_found=%d, nof_visited=%d", static_cast<int>(nodes.size()), nof_visited);
}

} // namespace mesh

---