MeshOutputFunctor.cpp

Go to the documentation of this file.

#include "MeshOutputFunctor.hpp"

#include "SplitVertices.hpp"
#include "moab/ParallelComm.hpp"
#include "RefinerTagManager.hpp"

#include <iostream>
#include <set>
#include <iterator>
#include <algorithm>

#undef MB_DEBUG

namespace moab {

MeshOutputFunctor::MeshOutputFunctor( RefinerTagManager* tag_mgr )
{
  this->mesh_in  = tag_mgr->get_input_mesh();
  this->mesh_out = tag_mgr->get_output_mesh();
  this->input_is_output = ( this->mesh_in == this->mesh_out );
  this->tag_manager = tag_mgr;
  this->destination_set = 0; // don't place output entities in a set by default.

  // When the input mesh and the output mesh are different, this map serves
  // as a dictionary from input vertices to output vertices.
  this->vertex_map = new SplitVertices<1>( this->tag_manager );

  // Hold information about newly-created vertices on subdivided edges and faces.
  this->split_vertices.resize( 4 );
  this->split_vertices[0] = 0; // Vertices (0-faces) cannot be split
  this->split_vertices[1] = new SplitVertices<1>( this->tag_manager );
  this->split_vertices[2] = new SplitVertices<2>( this->tag_manager );
  this->split_vertices[3] = new SplitVertices<3>( this->tag_manager );

  // Hold information about newly-created mesh entities (other than split vertices)
  // This is necessary in order for global IDs to be assigned consistently across processes.
  this->new_entities.resize( 5 );
  this->new_entities[0] = new EntitySource( 1, this->tag_manager );
  this->new_entities[1] = new EntitySource( 2, this->tag_manager );
  this->new_entities[2] = new EntitySource( 3, this->tag_manager );
  this->new_entities[3] = new EntitySource( 4, this->tag_manager );
  this->new_entities[4] = new EntitySource( 5, this->tag_manager );
}

MeshOutputFunctor::~MeshOutputFunctor()
{
  delete this->vertex_map;
  for ( int i = 1; i < 4; ++ i )
    delete this->split_vertices[i];
  for ( int i = 0; i < 5; ++ i )
    delete this->new_entities[i];
}

void MeshOutputFunctor::print_vert_crud(
    EntityHandle vout, int nvhash, EntityHandle* vhash, const double* vcoords, const void* /*vtags*/ )
{
  std::cout << "+ {";
  for ( int i = 0; i < nvhash; ++ i )
    std::cout << " " << vhash[i];
  std::cout << " } -> " << vout << " ";

  std::cout << "[ " << vcoords[0];
  for ( int i = 1; i < 6; ++i )
    std::cout << ", " << vcoords[i];
  std::cout << " ] ";

#if 0
  double* x = (double*)vtags;
  int* m = (int*)( (char*)vtags + 2 * sizeof( double ) );
  std::cout << "< " << x[0]
    << ", " << x[1];
  for ( int i = 0; i < 4; ++i )
    std::cout << ", " << m[i];
#endif // 0

  std::cout << " >\n";
  //std::cout << "##############################\n";
  //this->mesh_out->list_entities( 0, 1 );
  //std::cout << "##############################\n";
}

void MeshOutputFunctor::assign_global_ids( ParallelComm* comm )
{
  // First, we must gather the number of entities in each
  // partition (for all partitions, not just those resident locally).
  int lnparts = this->proc_partition_counts.size();
  std::vector<unsigned char> lpdefns;
  std::vector<int> lpsizes;
  lpdefns.resize( ProcessSet::SHARED_PROC_BYTES * lnparts );
  lpsizes.resize( lnparts );
#ifdef MB_DEBUG
  std::cout << "**** Partition Counts ****\n";
#endif // MB_DEBUG
  int i = 0;
  std::map<ProcessSet,int>::iterator it;
  for ( it = this->proc_partition_counts.begin(); it != this->proc_partition_counts.end(); ++ it, ++ i )
    {
    for ( int j = 0; j < ProcessSet::SHARED_PROC_BYTES; ++ j )
      lpdefns[ProcessSet::SHARED_PROC_BYTES * i + j] = it->first.data()[j];
    lpsizes[i] = it->second;
#ifdef MB_DEBUG
    std::cout << "Partition " << it->first << ": " << it->second << "\n";
#endif // MB_DEBUG
    }

  if ( ! comm )
    return;

  std::vector<int> nparts;
  std::vector<int> dparts;
  //unsigned long prank = comm->proc_config().proc_rank();
  unsigned long psize = comm->proc_config().proc_size();
  nparts.resize( psize );
  dparts.resize( psize + 1 );
  MPI_Allgather( &lnparts, 1, MPI_INT, &nparts[0], 1, MPI_INT, comm->proc_config().proc_comm() );
  //unsigned long ndefs = 0;
  for ( unsigned long rank = 1; rank <= psize; ++ rank )
    {
    dparts[rank] = nparts[rank - 1] + dparts[rank - 1];
#ifdef MB_DEBUG
    std::cout << "Proc " << rank << ": " << nparts[rank-1] << " partitions, offset: " << dparts[rank] << "\n";
#endif // MB_DEBUG
    }
  std::vector<unsigned char> part_defns;
  std::vector<int> part_sizes;
  part_defns.resize( ProcessSet::SHARED_PROC_BYTES * dparts[psize] );
  part_sizes.resize( dparts[psize] );
  MPI_Allgatherv(
    &lpsizes[0], lnparts, MPI_INT,
    &part_sizes[0], &nparts[0], &dparts[0], MPI_INT, comm->proc_config().proc_comm() );
  for ( unsigned long rank = 0; rank < psize; ++ rank )
    {
    nparts[rank] *= ProcessSet::SHARED_PROC_BYTES;
    dparts[rank] *= ProcessSet::SHARED_PROC_BYTES;
    }
  MPI_Allgatherv(
    &lpdefns[0], ProcessSet::SHARED_PROC_BYTES * lnparts, MPI_UNSIGNED_CHAR,
    &part_defns[0], &nparts[0], &dparts[0], MPI_UNSIGNED_CHAR, comm->proc_config().proc_comm() );

  // Now that we have the number of new entities in every partition, we
  // can deterministically assign the same GID to the same entity even
  // when shared across processors because we have an ordering that is
  // identical on all processes -- the vertex splits.
  for ( int j = 0; j < dparts[psize]; ++ j )
    {
    ProcessSet pset( &part_defns[ProcessSet::SHARED_PROC_BYTES * j] );
    std::map<ProcessSet,int>::iterator mit = this->proc_partition_counts.find( pset );
    if ( mit != this->proc_partition_counts.end() )
      {
#ifdef MB_DEBUG
      std::cout << "Partition " << pset << ( mit->second == part_sizes[j] ? " matches" : " broken" ) << ".\n";
#endif // MB_DEBUG
      }
    else
      {
      this->proc_partition_counts[pset] = part_sizes[j];
      }
    }
  
  std::map<ProcessSet,int> gids;
  std::map<ProcessSet,int>::iterator pcit;
  EntityHandle start_gid = 100; // FIXME: Get actual maximum GID across all processes and add 1
  for ( pcit = this->proc_partition_counts.begin(); pcit != this->proc_partition_counts.end(); ++ pcit )
    {
    gids[pcit->first] = start_gid;
    start_gid += pcit->second;
#ifdef MB_DEBUG
    std::cout << "Partition " << pcit->first << ": " << pcit->second << " # [" << gids[pcit->first] << "]\n";
#endif // MB_DEBUG
    }
  std::vector<SplitVerticesBase*>::iterator vit;
  vit = this->split_vertices.begin();
  ++ vit; // Skip split_vertices[0] since it's empty.
  ++ vit; // Skip split_vertices[1] since those entries already have global IDs... they exist in the input mesh.
  for ( /* skip */; vit != this->split_vertices.end(); ++ vit )
    {
    (*vit)->assign_global_ids( gids );
    }
  std::vector<EntitySource*>::iterator sit;
  for ( sit = this->new_entities.begin(); sit != this->new_entities.end(); ++ sit )
    {
    if ( *sit )
      (*sit)->assign_global_ids( gids );
    }
}

void MeshOutputFunctor::exchange_handles( ParallelComm*  )
{
}

void MeshOutputFunctor::assign_tags( EntityHandle vhandle, const void* vtags )
{
  if ( ! vhandle )
    return; // Ignore bad vertices

  int num_tags = this->tag_manager->get_number_of_vertex_tags();
  Tag tag_handle;
  int tag_offset;
  for ( int i = 0; i < num_tags; ++i )
    {
    this->tag_manager->get_output_vertex_tag( i, tag_handle, tag_offset );
    this->mesh_out->tag_set_data( tag_handle, &vhandle, 1, vtags );
    }
}

EntityHandle MeshOutputFunctor::map_vertex( EntityHandle vhash, const double* vcoords, const void* vtags )
{
  if ( this->input_is_output )
    { // Don't copy the original vertex!
#ifdef MB_DEBUG
    this->print_vert_crud( vhash, 1, &vhash, vcoords, vtags );
#endif // MB_DEBUG
    return vhash;
    }
  EntityHandle vertex_handle;
  bool newly_created = this->vertex_map->find_or_create(
    &vhash, vcoords, vertex_handle, this->proc_partition_counts, false );
  if ( newly_created )
    {
    std::vector<int> gid;
    this->assign_tags( vertex_handle, vtags );
    if ( this->tag_manager->get_input_gids( 1, &vhash, gid ) == MB_SUCCESS )
      {
      this->tag_manager->set_gid( vertex_handle, gid[0] );
      }
    }
  if ( ! vertex_handle )
    {
    std::cerr << "Could not insert vertex into new mesh!\n";
    }
#ifdef MB_DEBUG
  this->print_vert_crud( vertex_handle, 1, &vhash, vcoords, vtags );
  std::cout << "\nMap vert: " << vhash << " to: " << vertex_handle << "\n";
#endif // MB_DEBUG
  return vertex_handle;
}

EntityHandle MeshOutputFunctor::operator () ( int nvhash, EntityHandle* vhash, const double* vcoords, const void* vtags )
{
  EntityHandle vertex_handle;
  if ( nvhash < 4 )
    {
    bool newly_created = this->split_vertices[nvhash]->find_or_create(
      vhash, vcoords, vertex_handle, this->proc_partition_counts, true );
    if ( newly_created )
      {
      this->assign_tags( vertex_handle, vtags );
      }
    if ( ! vertex_handle )
      {
      std::cerr << "Could not insert mid-edge vertex!\n";
      }
#ifdef MB_DEBUG
    std::cout << "(-" << nvhash << "-) ";
    this->print_vert_crud( vertex_handle, nvhash, vhash, vcoords, vtags );
#endif // MB_DEBUG
    }
  else
    {
    vertex_handle = 0;
    std::cerr << "Not handling splits on faces with " << nvhash << " corners yet.\n";
    }
  return vertex_handle;
}

void MeshOutputFunctor::operator () ( EntityHandle h )
{
#ifdef MB_DEBUG
  std::cout << h << " ";
#endif // MB_DEBUG
  if ( ! this->input_is_output )
    {
    // FIXME: Copy to output mesh
    }
  this->elem_vert.push_back( h );
}

void MeshOutputFunctor::operator () ( EntityType etyp )
{
  EntityHandle elem_handle;
  int nconn = this->elem_vert.size();
  bool newly_created = this->new_entities[nconn]->create_element(
    etyp, nconn, &this->elem_vert[0], elem_handle, this->proc_partition_counts );
  if ( newly_created )
    {
#ifdef MB_DEBUG
    std::cout << " *** ";
#endif // MB_DEBUG
    // FIXME: Handle tag assignment for elements as well as vertices
    this->tag_manager->assign_element_tags( elem_handle );
    }
#ifdef MB_DEBUG
  std::cout << "---------> " << elem_handle << " ( " << etyp << " )\n\n";
#endif // MB_DEBUG
  this->elem_vert.clear();
}

} // namespace moab