src/tetmesh/paralleltetmesh.cpp

Go to the documentation of this file.

//
// C++ Interface: parallel_tetmesh
//
// Description: 
//
//
// Author: Dag Evensberget <deg.evensberget@psi.ch>, (C) 2005
//
// Copyright: See COPYING file that comes with this distribution
//
//

#include <set>
#include <stdexcept>
#include "myrlog.h"
#include "paralleltetmesh.h"

using namespace std;

namespace mesh {

ParallelTetMesh::ParallelTetMesh(MPI_Comm comm)
    : TetMesh(), _num_global_edges(0) 
{
    _comm = comm;   
    MPI_Comm_size(_comm, &_mpi_size);
    MPI_Comm_rank(_comm, &_mpi_rank);

    _tet_counter = 0;
    _point_counter = 0;
}

ParallelTetMesh::~ParallelTetMesh() {
}
        
void ParallelTetMesh::initPoint(int nofPoint) {
    _point_gids.clear();
    _point_counter = 0;
    TetMesh::initPoint(nofPoint);
 }   
 
void ParallelTetMesh::initTet(int nofTet) {
    _tet_counter = 0;    
    TetMesh::initTet(nofTet);
    
/*    for(int i = 0; i < _point_counter; i ++) 
         cerr<<"G->L "<<_point_global_id[i]<<"->"<<_point_local_id[_point_global_id[i]]<<endl;*/
   
}

void ParallelTetMesh::insertPoint(int global_id, double x, double y, double z) {
    pair<int,int> p(global_id, _point_counter); 
    _point_local_id.insert(p);
    //cerr << "Proc " << _mpi_rank << " inserting point g->l " << global_id << "->" << _point_counter << endl;    
    _point_gids.push_back(global_id);    
    TetMesh::insertPoint(_point_counter, x, y, z);
    _point_counter += 1;
}

  void ParallelTetMesh::insertTet(id_t global_id, id_t id0, id_t id1, id_t id2, id_t id3, id_t material){
    pair<int,int> p(global_id, _tet_counter);
    //Insert with local numbers
    TetMesh::insertTet(_tet_counter, _point_local_id[id0], _point_local_id[id1], _point_local_id[id2], _point_local_id[id3], material);
    _tet_counter += 1;
}
   
id_t ParallelTetMesh::lookupFace(id_t node0, id_t node1, id_t node2) {
//        cerr <<"Proc " << _mpi_rank << " looking up g->l "<< node0<<"->" << _point_local_id[node0]<<" "<< node1<<"->"<< //_point_local_id[node1]<<" "<< node2<<"->"<< _point_local_id[node2]<<endl;
    
    map<int,int>::iterator iter;
    iter = _point_local_id.find(node0);
    if (iter == _point_local_id.end())
        throw std::runtime_error("Accessed non-existing face (in ParallelTetMesh).");
    iter = _point_local_id.find(node1);
    if (iter == _point_local_id.end())
        throw std::runtime_error("Accessed non-existing face (in ParallelTetMesh).");
    iter = _point_local_id.find(node2);
    if (iter == _point_local_id.end())
        throw std::runtime_error("Accessed non-existing face (in ParallelTetMesh).");
        
    id_t face = TetMesh::lookupFace(_point_local_id[node0], _point_local_id[node1], _point_local_id[node2]);
    return face;
}

void ParallelTetMesh::fix_edge_orientation() {
    for (tet_iterator iter = tet_begin(); iter != tet_end(); ++ iter) {
        for (int j = 0; j < 6; j ++) {
            // get local edge j
            Edge* edge = (*iter).get_edge(j);
            rAssert(edge->get_start_id() < edge->get_end_id());
            
            // test global orientation
            if (point_gid(edge->get_start_id()) >= point_gid(edge->get_end_id())) {
                rWarningAll("Global edge orientation error.");
                rAssert(false);
                // FIXME: Should swap edge orientation here
            }
        }
    }
}

void ParallelTetMesh::finalize_mesh(int nofSym) {
    synchronize();
    TetMesh::finalize_mesh(nofSym);
}

void ParallelTetMesh::synchronize() {

    // Assign ownership to nodes (the node ownership info
    set_point_ownership();

    // Assign global ids to faces
    set_face_gids();
    
    // Assign global ids to edges
    set_edge_gids();
    
    // Fix edge orientation (must be done consistently across all processors)
    fix_edge_orientation();
}

utility::triple<id_t, id_t, id_t> ParallelTetMesh::get_face_corner_gids(const Face& face) const {
    id_t c0 = point_gid(face.get_corner_id(0));
    id_t c1 = point_gid(face.get_corner_id(1));
    id_t c2 = point_gid(face.get_corner_id(2));
    if (c0 > c1)
        std::swap<id_t>(c0, c1);
    if (c1 > c2)
        std::swap<id_t>(c1, c2);
    if (c0 > c1)
        std::swap<id_t>(c0, c1);                    
    return utility::triple<id_t, id_t, id_t>(c0, c1, c2);
}

std::pair<id_t, id_t> ParallelTetMesh::get_edge_corner_gids(const Edge& edge) const {
    id_t c0 = point_gid(edge.get_start_id());
    id_t c1 = point_gid(edge.get_end_id());
    if (c0 > c1)
        std::swap<id_t>(c0, c1);
    return std::pair<id_t, id_t>(c0, c1);
}

void ParallelTetMesh::set_face_gids() {
    // Find interprocessor boundary faces
    typedef utility::triple<id_t, id_t, id_t> ipb_face_t;
    vector<ipb_face_t> bfaces_global_node_ids;
    for (face_iterator iter = face_begin(); iter != face_end(); ++iter) {
        id_t tet0, tet1;
        (*iter).get_tet_ids(&tet0, &tet1);
        if (tet1 == ID_NONE) {
            if (!(*iter).on_boundary() && (*iter).on_symmetry_planes() == 0) {
                //Face is on iterprocessor boundary.
                bfaces_global_node_ids.push_back(get_face_corner_gids(*iter));
            }
        }
    }
    // Assign global ids for ipb faces     
    vector<id_t> ipb_face_gids;
    vector<int> ipb_faces_counts, ipb_faces_displs;
    int num_global_ipb_faces = parallel::assign_gids(bfaces_global_node_ids, 
                                                     ipb_face_gids,
                                                     ipb_faces_counts,
                                                     ipb_faces_displs,
                                                     _comm);
    
    //
    // Assign gids for all faces (in order of local face ids)
    //
    int num_my_nonipb_faces = get_nof_faces() - bfaces_global_node_ids.size();
    std::vector<int> counts, displs;
    int num_global_nonipb_faces = \
        parallel::set_counts_displs(num_my_nonipb_faces, counts, displs, _comm);
    int non_ipb_face_counter = 0;
    for (face_iterator iter = face_begin(); iter != face_end(); ++ iter) {
        id_t tet0, tet1;
        (*iter).get_tet_ids(&tet0, &tet1);
        if ((tet1 == ID_NONE) &&
            (!(*iter).on_boundary() && (*iter).on_symmetry_planes() == 0)) {
            //
            // Face is on iterprocessor boundary.
            //
            pair<vector<ipb_face_t>::iterator, vector<ipb_face_t>::iterator> result = \
                equal_range(bfaces_global_node_ids.begin(), 
                            bfaces_global_node_ids.end(), 
                            get_face_corner_gids(*iter));
            assert(result.second - result.first == 1); // Want to find exactly one element.

            // Assign corresponding ipb gid            
            _face_gids.push_back(ipb_face_gids[result.first - bfaces_global_node_ids.begin()]);
        } else {
            //
            // Inner face or face on geom boundary or symmetry plane
            //
            _face_gids.push_back(non_ipb_face_counter + num_global_ipb_faces + displs[_mpi_rank]);
            non_ipb_face_counter += 1;
        }
    }
    assert(non_ipb_face_counter == num_my_nonipb_faces);
    assert(_face_gids.size() == get_nof_faces());
    _num_global_faces = num_global_ipb_faces + num_global_nonipb_faces;
    rInfoAll("Number of faces: (owned: %d, local: %d, global: %d)", 
             ipb_faces_counts[_mpi_rank]+num_my_nonipb_faces, get_nof_faces(), _num_global_faces);
    rInfoAll("Owned faces: ipb:[%d:%d), non-ipb:[%d:%d)",
             ipb_faces_displs[_mpi_rank], 
             ipb_faces_displs[_mpi_rank] + ipb_faces_counts[_mpi_rank],
             num_global_ipb_faces + displs[_mpi_rank], 
             num_global_ipb_faces + displs[_mpi_rank] + counts[_mpi_rank]);

    _my_faces[0] = ipb_faces_displs[_mpi_rank];
    _my_faces[1] = ipb_faces_displs[_mpi_rank] + ipb_faces_counts[_mpi_rank];
    _my_faces[2] = num_global_ipb_faces + displs[_mpi_rank];
    _my_faces[3] = num_global_ipb_faces + displs[_mpi_rank] + counts[_mpi_rank];
    
    // debug
#if 0
    vector<id_t> eliminated_gids;
    for (face_iterator iter = face_begin(); iter != face_end(); ++ iter) {
        if ((*iter).on_boundary())
            eliminated_gids.push_back(face_gid((*iter).get_id()));
    }
    std::map<id_t,id_t> mgid_map;
    id_t num_global_mgids;
    parallel::compute_mgid_map(num_global_faces,
                               _face_gids,
                               eliminated_gids,
                               &num_global_mgids,
                               mgid_map,
                               _comm);
    rInfoAll("num_global_mgids=%d", num_global_mgids);
#endif
}

void ParallelTetMesh::set_edge_gids() {
    // Find interprocessor boundary edge
    typedef std::pair<id_t, id_t> ipb_edge_t;    
    vector<ipb_edge_t> bedges_global_node_ids;
    set<id_t> edges_ipb;
    for (face_iterator iter = face_begin(); iter != face_end(); ++iter) {
        id_t tet0, tet1;
        (*iter).get_tet_ids(&tet0, &tet1);
        if (tet1 == ID_NONE) {
            if (!(*iter).on_boundary() && (*iter).on_symmetry_planes() == 0) {
                /*
                 * Face is on iterprocessor boundary. Check if edge has already been added.
                 * If not, add to bedges_global_node_ids.
                 */
                for (int i = 0; i < 3; i++) {                    
                    Edge* e = (*iter).get_edge(i);               
                    if (edges_ipb.find(e->get_id()) == edges_ipb.end()) {
                        bedges_global_node_ids.push_back(get_edge_corner_gids(*e));
#if 1
                    ipb_edge_t t = get_edge_corner_gids(*e);
                    if (t.first == 376 && t.second == 589) {
                        rWarningAll("edge (376, 589) on ipb");
                    }
#endif
                        edges_ipb.insert(e->get_id());
                    } else {
                        // Edge has already been added.
                    }
                }                
            }
        }
    }
    
    // For each ipb edge decide whether this processor should claim it and store in boolean array
    bool* edges_claim = new bool[bedges_global_node_ids.size()];
    for (size_t i = 0; i < bedges_global_node_ids.size(); ++ i)
    {
        // Only claim edge if at least one of its corners are owned by this processor.
        // This guarantees a ML-compatible matrix distribution.
        ipb_edge_t* it = &(bedges_global_node_ids[i]);
        edges_claim[i] = is_owned_point((*it).first) || is_owned_point((*it).second);
    }
        
    rDebugAll("bedges_global_node_ids.size() = %d", 
              static_cast<int>(bedges_global_node_ids.size()));
    
    // Assign global ids for ipb edges.
    vector<id_t> ipb_edge_gids;
    vector<int> ipb_edges_counts, ipb_edges_displs;
    
    /* FIXME: assign_gids2 does not work in all cases. If we are unlucky a processor gets to own
     *        both points of an edge but does not actually store that edge locally. In this 
     *        assign_gids is unable to find a proper owner for that edge and terminates with
     *        a failed assertion.
     * 
     *        It is not clear whether a proper distribution of edges and points satisfying the
     *        ML constraint exists in all cases.
     * 
     *        One solution would be to allow a processor to own an edge even if it does not
     *        store it locally. Some portions of the mesh code may need to be adapted to
     *        implement this.
     * 
     *        For the moment the old assign_gids is used to determine ownership of the edges.
     *        The code will thus not generate ML compatible distributions. However it will 
     *        generate a distribution which should work in all cases if Maxwell ML is not used.
     */
    //int num_global_ipb_edges = parallel::assign_gids2(edges_claim, bedges_global_node_ids, ipb_edge_gids, ipb_edges_counts, ipb_edges_displs, _comm);
    int num_global_ipb_edges = parallel::assign_gids(bedges_global_node_ids, ipb_edge_gids, ipb_edges_counts, ipb_edges_displs, _comm);
    
    delete edges_claim;
    
    // debugging
#if 0
    for (int i = 0; i < bedges_global_node_ids.size(); ++ i)
        rDebugAll("ipb_edge_gids[%d] = %d", i, ipb_edge_gids[i]);
#endif
    
    //
    // Assign global ID's for all edges. In local edge order.
    // 
    int num_my_nonipb_edges = get_nof_edges() - bedges_global_node_ids.size();
    std::vector<int> counts, displs;
    int num_global_nonipb_edges = \
        parallel::set_counts_displs(num_my_nonipb_edges, counts, displs, _comm);
    int non_ipb_edge_counter = 0;
    
    for (edge_iterator it = edge_begin(); it != edge_end(); ++it) {
        if (edges_ipb.find((*it).get_id()) != edges_ipb.end()) {
            //
            // Edge is on interprocessor boundary            
            // 
            pair<vector<ipb_edge_t>::iterator, vector<ipb_edge_t>::iterator> result = \
                equal_range(bedges_global_node_ids.begin(), bedges_global_node_ids.end(), get_edge_corner_gids(*it));
            assert (result.second - result.first == 1);
            // Assign corresponding ipb gid.
            _edge_gids.push_back(ipb_edge_gids[result.first - bedges_global_node_ids.begin()]);
        } else {
            //
            // Inner edge or bondary edge on geometrical boundary or symmetry plane.
            //
            _edge_gids.push_back(non_ipb_edge_counter + num_global_ipb_edges + displs[_mpi_rank]);
            non_ipb_edge_counter++;
        }
    }
    
    assert(non_ipb_edge_counter == num_my_nonipb_edges);
    assert(_edge_gids.size() == get_nof_edges());
    _num_global_edges = num_global_ipb_edges + num_global_nonipb_edges;
    rInfoAll("Number of edges: (owned: %d, local: %d, global: %d)", 
             ipb_edges_counts[_mpi_rank]+num_my_nonipb_edges, get_nof_edges(), _num_global_edges);
    rInfoAll("Owned edges: ipb:[%d:%d), non-ipb:[%d:%d)",
             ipb_edges_displs[_mpi_rank], 
             ipb_edges_displs[_mpi_rank] + ipb_edges_counts[_mpi_rank],
             num_global_ipb_edges + displs[_mpi_rank], 
             num_global_ipb_edges + displs[_mpi_rank] + counts[_mpi_rank]);
    _my_edges[0] = ipb_edges_displs[_mpi_rank];
    _my_edges[1] = ipb_edges_displs[_mpi_rank] + ipb_edges_counts[_mpi_rank];
    _my_edges[2] = num_global_ipb_edges + displs[_mpi_rank];
    _my_edges[3] = num_global_ipb_edges + displs[_mpi_rank] + counts[_mpi_rank];
   
    
#if 0             
    // debug
    vector<id_t> eliminated_gids;
    for (edge_iterator iter = edge_begin(); iter != edge_end(); ++iter) {
        if ((*iter).on_boundary())
            eliminated_gids.push_back(edge_gid((*iter).get_id()));
    }
    std::map<id_t,id_t> mgid_map;
    id_t num_global_mgids;
    parallel::compute_mgid_map(_num_global_edges,
                               _edge_gids,
                               eliminated_gids,
                               &num_global_mgids,
                               mgid_map,
                               _comm);
    rInfoAll("num_global_mgids=%d", num_global_mgids);
#endif    
}

void ParallelTetMesh::set_point_ownership() {
    int mpi_rank;
    MPI_Comm_rank(_comm, &mpi_rank);

    //
    // Determine ipb points
    //
    set<id_t> ipb_point_gids_set;
    // Loop over ipb faces
    for (face_iterator iter = face_begin(); iter != face_end(); ++ iter) {
        id_t tet0, tet1;
        (*iter).get_tet_ids(&tet0, &tet1);
        if ((tet1 == ID_NONE) &&
            (!(*iter).on_boundary() && (*iter).on_symmetry_planes() == 0)) {
            //
            // Face is on iterprocessor boundary.
            //
            id_t c0 = point_gid((*iter).get_corner_id(0));
            id_t c1 = point_gid((*iter).get_corner_id(1));
            id_t c2 = point_gid((*iter).get_corner_id(2));
            ipb_point_gids_set.insert(c0);
            ipb_point_gids_set.insert(c1);
            ipb_point_gids_set.insert(c2);
        }
    }    
    vector<id_t> ipb_point_gids;
    // Store unique point gids in ascending order
    for (set<id_t>::iterator it = ipb_point_gids_set.begin();
         it != ipb_point_gids_set.end();
         ++ it)
    {
        ipb_point_gids.push_back(*it);
    }
    ipb_point_gids_set.clear();
    
    // Set ownership of ipc points
    vector<int> ipb_point_owners;
    int num_ipb_points_owned = parallel::set_ownership(ipb_point_gids, ipb_point_owners, _comm);
    
    // Set ownership of all locally stored points
    int num_points_owned = num_ipb_points_owned;    
    _point_owned_by_me.resize(get_nof_points());
    for (point_iterator iter = point_begin(); iter != point_end(); ++ iter) {
        id_t lid = (*iter).get_id();
        id_t gid = point_gid(lid);
        vector<id_t>::iterator found = \
            lower_bound(ipb_point_gids.begin(), ipb_point_gids.end(), gid);
        if (found == ipb_point_gids.end() || *found != gid) {
            // Not on ipb
            _point_owned_by_me[lid] = true;
            num_points_owned += 1;
        } else {
            // On ipb
            if (ipb_point_owners[found - ipb_point_gids.begin()] == mpi_rank)
                _point_owned_by_me[lid] = true;
            else
                _point_owned_by_me[lid] = false;
        }
    }
    
    //Set nof points globally    
    MPI_Allreduce (&num_points_owned, &_num_global_points, 1, MPI_INT, MPI_SUM, _comm);
    
    rDebugAll("Owned points: ipb:%d, total:%d.", num_ipb_points_owned, num_points_owned);
    rDebug("Number of points globally: %d", _num_global_points);
    
    // debugging
#if 0
    for (point_iterator iter = point_begin(); iter != point_end(); ++ iter) {
        id_t lid = (*iter).get_id();
        if (_point_owned_by_me[lid])
            rDebugAll("own point %d", point_gid(lid));
    }
#endif    
}

id_t ParallelTetMesh::face_gid(id_t lid) const {
    return _face_gids[lid];
}

id_t ParallelTetMesh::edge_gid(id_t lid) const {
    return _edge_gids[lid];
}

id_t ParallelTetMesh::point_gid(id_t lid) const {
    return _point_gids[lid];
}

id_t ParallelTetMesh::get_num_global_tets()
{
        int i,j,k; 
        id_t ntet; 
        id_t gntet; 
        ntet=get_nof_tets();
        MPI_Allreduce(&ntet,&gntet,1,MPI_UNSIGNED,MPI_SUM,_comm);

        /* go back */
        return(gntet);
}




bool ParallelTetMesh::is_owned_point(id_t gid) const {
    map<int,int>::const_iterator it = _point_local_id.find(gid);
    if (it == _point_local_id.end())
        return false;
    else
        return _point_owned_by_me[(*it).second];
}

bool ParallelTetMesh::is_owned_edge(id_t gid) const {
    return (_my_edges[0] <= gid && gid < _my_edges[1]) \
        || (_my_edges[2] <= gid && gid < _my_edges[3]);
}

bool ParallelTetMesh::is_owned_face(id_t gid) const {
    return (_my_faces[0] <= gid && gid < _my_faces[1]) \
        || (_my_faces[2] <= gid && gid < _my_faces[3]);
}

MPI_Comm ParallelTetMesh::get_comm() {
    return _comm;
}

} // namespace mesh

---