//#include <mpi.h>
#include <moab/Core.hpp>
#include <moab/Interface.hpp>
#include <moab/ParallelComm.hpp>
#include <moab/MeshTopoUtil.hpp>
#include <moab/CN.hpp>
#include <MBZoltan.hpp>
#include <MBTypes.h>
#include <MBTagConventions.hpp>
#include <MBParallelConventions.h>
#include <sstream>
#include <fstream>
#include <vector>
#include <metis.h>
#include <parmetis.h>

//#include <petsc.h>

//#define USE_EXCHANGE_OWNED_MESHES
//#define USE_SENDRECV
//#define DEBUG_OUTPUT
//#define STORE_REMOTE_HANDLES

#define CHKERR(...) {if(result) {std::string str; imoab.get_last_error(str); std::cout << pcomm.rank() << " line: " << __LINE__ << " reason: " << str << std::endl; printf(__VA_ARGS__); MPI_Abort(MPI_COMM_WORLD,__LINE__);}}
#define FAIL {std::cout << icomm.rank() << ":" << __FILE__ << ":" << __FUNCTION__ << ":" << __LINE__ << std::endl; return MB_FAILURE;}
#define RFAIL {std::cout << icomm.rank() << ":" << __FILE__ << ":" << __FUNCTION__ << ":" << __LINE__ << std::endl; return result;}

void start_serial_output(MPI_Comm comm)
{
	int rank;
	MPI_Barrier(comm);
	MPI_Comm_rank(comm,&rank);
	for(int i = 0; i < rank; i++) MPI_Barrier(comm);
}

void end_serial_output(MPI_Comm comm)
{
	int rank,size;
	MPI_Comm_rank(comm,&rank);
	MPI_Comm_size(comm,&size);
	for(int i = rank; i < size; i++) MPI_Barrier(comm);
}

int PIMisc_delete_ghost(moab::Interface & imoab,EntityHandle set_handle, moab::ParallelComm & icomm,int elem_dim)
{
	moab::Range ents, keep, erase;
	int result,dims , i;
	unsigned char pstatus;
	ents.clear();
	result = imoab.get_entities_by_dimension(set_handle,elem_dim,ents);
	if( result != MB_SUCCESS) RFAIL;
	for(moab::Range::iterator ite = ents.begin(); ite != ents.end(); ite++)
	{
		result = imoab.tag_get_data(icomm.pstatus_tag(),&*ite,1,&pstatus);
		if( result != MB_SUCCESS) RFAIL;
		if( pstatus & (PSTATUS_GHOST) )
			imoab.delete_entities(&*ite,1);
	}
	result = imoab.get_dimension(dims);
	if( result != MB_SUCCESS ) RFAIL;
	ents.clear();
	result = imoab.get_entities_by_dimension(set_handle,elem_dim,ents);
	if( result != MB_SUCCESS) RFAIL;
	for(i = dims; i >= 0; i--)
	{
		keep.clear();
		erase.clear();
		if( i != elem_dim )
		{
			result = imoab.get_entities_by_dimension(set_handle,i,erase);
			if( result != MB_SUCCESS) RFAIL;
			result = imoab.get_adjacencies(ents,i,0,keep,moab::Interface::UNION);
			if( result != MB_SUCCESS) RFAIL;
			erase -= keep;
			imoab.delete_entities(erase);
		}
	}
	return MB_SUCCESS;	
}


int PIFile_write_pvtk(moab::Interface & imoab, EntityHandle set_handle, moab::ParallelComm & icomm, std::string filename)
{
	//TASK don't delete ghost entities, create set of owned entities and save this set
	int result;
	PIMisc_delete_ghost(imoab,imoab.get_root_set(),icomm,3);
	std::stringstream save_file;
	save_file << filename << icomm.rank() << ".vtk";
	result = imoab.write_file(save_file.str().c_str());
	if( result != MB_SUCCESS ) RFAIL;
	if( icomm.rank() == 0 )
	{
		std::stringstream save_pvtk_file;
		std::fstream pvtk;
		save_pvtk_file << filename << ".pvtk";
		pvtk.open(save_pvtk_file.str().c_str(),std::fstream::out);
		pvtk << "<File version=\"pvtk-1.0\" dataType=\"vtkUnstructuredGrid\" numberOfPieces=\"" 
			 << icomm.size() << "\">" << std::endl;
		for(unsigned int i = 0; i < icomm.size(); i++)
			pvtk << "<Piece fileName=\"" << filename << i << ".vtk\"/>" << std::endl;
		pvtk << "</File>" << std::endl;
		pvtk.close();
	}
	return MB_SUCCESS;
}

/*
int PIFile_read_pvtk(moab::Interface & imoab, EntityHandle set_handle, moab::ParallelComm & icomm, std::string filename)
{
	int result;
	if( result != MB_SUCCESS ) RFAIL;
	if( icomm.rank() == 0 )
	{
		std::stringstream save_pvtk_file;
		std::fstream pvtk;
		save_pvtk_file << filename << ".pvtk";
		pvtk.open(save_pvtk_file.str().c_str(),std::fstream::out);
		pvtk << "<File version=\"pvtk-1.0\" dataType=\"vtkUnstructuredGrid\" numberOfPieces=\"" 
			 << icomm.size() << "\">" << std::endl;
		for(unsigned int i = 0; i < icomm.size(); i++)
			pvtk << "<Piece fileName=\"" << filename << i << ".vtk\"/>" << std::endl;
		pvtk << "</File>" << std::endl;
		pvtk.close();
	}
	return MB_SUCCESS;
}
*/


int PIFile_load_distributed(moab::Interface & imoab, EntityHandle set_handle, moab::ParallelComm & icomm, char * filename)
{
	int result;
	std::stringstream options;
	options << "PARALLEL=READ_PART;";
	options << "PARTITION=PARALLEL_PARTITION;";
	options << "PARTITION_BY_RANK;";
	options << "PARALLEL_COMM=" << icomm.get_id() << ";";
	result = imoab.load_file(filename,set_handle == imoab.get_root_set() ? NULL : &set_handle,options.str().c_str());
	if( result != MB_SUCCESS ) RFAIL;
	return result;
}

//loads file on 1 processor
int PIFile_load_serial_rank(moab::Interface & imoab,EntityHandle set_handle, moab::ParallelComm & icomm, char * filename, unsigned int rank)
{
	int temp = 0, dim,result;
	if( icomm.rank() == rank )
	{
		int dims;
		//int i,j;
		//moab::Tag global_id;
		moab::Range elems;
		std::vector<int> ids;
		//load file
		{
			result = imoab.load_file(filename,set_handle == imoab.get_root_set() ? NULL : &set_handle);
			if( result != MB_SUCCESS ) RFAIL;
		}
		//get number of dimensions
		{
			result = imoab.get_dimension(dims);
			if( result != MB_SUCCESS ) RFAIL;
		}
		temp = dim;
	}
	//share number of dimensions
	{
		result = MPI_Allreduce(&temp,&dim,1,MPI_INTEGER,MPI_MIN,icomm.comm());
		if( result != MPI_SUCCESS ) FAIL;
		if( !dim ) dim = 3;
		result = imoab.set_dimension(dim);
		if( result != MB_SUCCESS ) RFAIL;
	}
	//delete any already existing PARALLEL_PARTITION set
	{
		int i;
		moab::Range parallel_sets,elems;
		moab::Tag parallel_partition = icomm.part_tag();
		result = imoab.get_entities_by_type_and_tag(set_handle,moab::MBENTITYSET,&parallel_partition,NULL,1,parallel_sets,moab::Interface::UNION);
		if( result != MB_SUCCESS ) RFAIL;
		if( !parallel_sets.empty() )
		{
			result = imoab.clear_meshset(parallel_sets);
			if( result != MB_SUCCESS ) RFAIL;
		}
		if( parallel_sets.size() == 0 )
		{
			EntityHandle new_set;
			result = imoab.create_meshset(moab::MESHSET_SET,new_set);
			if( result != MB_SUCCESS ) RFAIL;
			parallel_sets.insert(new_set);
		}
		else if( parallel_sets.size() != 1 )
		{
			int num_del = parallel_sets.size()-1;
			for(i = 0; i < num_del; i++)
			{
				EntityHandle old_set = parallel_sets.pop_back();
				result = imoab.delete_entities(&old_set,1);
				if( result != MB_SUCCESS ) RFAIL;
			}
		}
		result = imoab.tag_set_data(parallel_partition,parallel_sets,&rank);
		if( result != MB_SUCCESS ) RFAIL;
		for( i = dim; i >= dim; i--)
		{
			result = imoab.get_entities_by_dimension(set_handle,i,elems,1);
			if( result != MB_SUCCESS ) RFAIL;
			result = imoab.add_entities(parallel_sets[0],elems);
			if( result != MB_SUCCESS ) RFAIL;
			elems.clear();
		}
		result = imoab.get_entities_by_type_and_tag(set_handle,MBENTITYSET,&parallel_partition,NULL,1,icomm.partition_sets());
		if( result != MB_SUCCESS ) RFAIL;
	}
	return MB_SUCCESS;
}

//loads file and cut between processors
int PIFile_load_serial(moab::Interface & imoab,EntityHandle set_handle, moab::ParallelComm & icomm, char * filename)
{
	moab::Tag parallel_partition;
	moab::Tag global_id;
	int dims,result = MB_SUCCESS,i;
	unsigned int local_start,local_end,local_size,rank = icomm.rank();
	unsigned int j;
	moab::Range erase[4], keep[4], parallel_sets, elems;
	//load mesh on every processor
	{
		result = imoab.load_file(filename,set_handle == imoab.get_root_set() ? NULL : &set_handle);
		if( result != MB_SUCCESS ) RFAIL;
	}
	//get number of mesh dimensions
	{
		result = imoab.get_dimension(dims);
		if( result != MB_SUCCESS ) RFAIL;
	}
	//assign id
	{
		moab::Range ents;
		std::vector<int> ids;
		result = imoab.tag_get_handle(PARALLEL_GID_TAG_NAME,1,moab::MB_TYPE_INTEGER,global_id);
		if( result != MB_SUCCESS ) RFAIL;
		for(i = dims; i >= 0; i--)
		{
			result = imoab.get_entities_by_dimension(set_handle,i,ents);
			if( result != MB_SUCCESS ) RFAIL;
			ids.resize(ents.size());
			for(j = 0; j < ents.size(); j++) ids[j] = j;
			result = imoab.tag_set_data(global_id,ents,&ids[0]);
			if( result != MB_SUCCESS ) RFAIL;
		}
	}
	//collect local data
	{
		result = imoab.get_entities_by_dimension(set_handle,dims,erase[dims]);
		if( result != MB_SUCCESS ) RFAIL;
		local_size = erase[dims].size()/icomm.size();
		if( local_size == 0 ) FAIL;
		local_start = local_size*icomm.rank();
		local_end = local_size*(icomm.rank()+1);
		if( icomm.rank() == icomm.size()-1 ) local_end = erase[dims].size();
		for(i = local_start; i < (int)local_end; i++)
			keep[dims].insert(erase[dims][i]);
		erase[dims] -= keep[dims];
		for(i = dims-1; i >= 0; i--)
		{
			result = imoab.get_adjacencies(keep[i+1],i,1,keep[i],moab::Interface::UNION);
			if( result != MB_SUCCESS ) RFAIL;
			result = imoab.get_entities_by_dimension(set_handle,i,erase[i]);
			if( result != MB_SUCCESS ) RFAIL;
			erase[i] -= keep[i];
		}
	}
	//erase non-local data
	{
		for(i = dims; i >= 0; i--)
		{
			result = imoab.delete_entities(erase[i]);
			if( result != MB_SUCCESS ) RFAIL;
		}
	}
	//assign PARALLEL_PARTITION tag for set of local entities
	//delete any already existing PARALLEL_PARTITION set
	{
		//result = imoab.tag_get_handle(PARALLEL_PARTITION_TAG_NAME,1,moab::MB_TYPE_INTEGER,parallel_partition,moab::MB_TAG_SPARSE|moab::MB_TAG_CREAT);
		//if( result != MB_SUCCESS ) RFAIL;
		parallel_partition = icomm.part_tag();
		result = imoab.get_entities_by_type_and_tag(set_handle,moab::MBENTITYSET,&parallel_partition,NULL,1,parallel_sets,moab::Interface::UNION);
		if( result != MB_SUCCESS ) RFAIL;
		if( !parallel_sets.empty() )
		{
			result = imoab.clear_meshset(parallel_sets);
			if( result != MB_SUCCESS ) RFAIL;
		}
		if( parallel_sets.size() == 0 )
		{
			EntityHandle new_set;
			result = imoab.create_meshset(moab::MESHSET_SET,new_set);
			if( result != MB_SUCCESS ) RFAIL;
			parallel_sets.insert(new_set);
		}
		else if( parallel_sets.size() != 1 )
		{
			int num_del = parallel_sets.size()-1;
			for(i = 0; i < num_del; i++)
			{
				EntityHandle old_set = parallel_sets.pop_back();
				result = imoab.delete_entities(&old_set,1);
				if( result != MB_SUCCESS ) RFAIL;
			}
		}
		result = imoab.tag_set_data(parallel_partition,parallel_sets,&rank);
		if( result != MB_SUCCESS ) RFAIL;
		for( i = dims; i >= dims; i--)
		{
			result = imoab.get_entities_by_dimension(set_handle,i,elems,1);
			if( result != MB_SUCCESS ) RFAIL;
			result = imoab.add_entities(parallel_sets[0],elems);
			if( result != MB_SUCCESS ) RFAIL;
			elems.clear();
		}
		result = imoab.get_entities_by_type_and_tag(set_handle,MBENTITYSET,&parallel_partition,NULL,1,icomm.partition_sets());
		if( result != MB_SUCCESS ) RFAIL;
	}	
	return MB_SUCCESS;
}



int PIMisc_split_entities(moab::Interface & imoab, EntityHandle set_handle, moab::ParallelComm & icomm, moab::Range & local, moab::Range & ghost, int elem_dim)
{
	int result,owner,rank;
	moab::Range::iterator obj;
	result = imoab.get_entities_by_dimension(set_handle,elem_dim,local);
	if( result != MB_SUCCESS ) RFAIL;
	obj = local.begin();
	rank = icomm.rank();
	while(obj != local.end())
	{
		result = icomm.get_owner(*obj,owner);
		if( result != MB_SUCCESS ) RFAIL;
		if( owner == rank ) 
			obj++;
		else
		{
			ghost.insert(*obj);
			obj = local.erase(obj);
		}
	}
	return MB_SUCCESS;
}

#define PIDIST_PART_WEIGHT_TAG_NAME "PIDIST_PART_WGT"
#define PIDIST_PART_WEIGHT_TYPE moab::MB_TYPE_DOUBLE
#define PIDIST_PART_WEIGHT_TOLERANCE_TAG_NAME "PIDIST_WEIGHT_TOLERANCE"
#define PIDIST_PART_WEIGHT_TOLERANCE_TYPE moab::MB_TYPE_DOUBLE
#define PIDIST_VERTEX_WEIGHT_TAG_NAME "PIDIST_VERTEX_WEIGHT"
#define PIDIST_VERTEX_WEIGHT_TYPE moab::MB_TYPE_INTEGER
#define PIDIST_ADJACENCIES_WEIGHT_TAG_NAME "PIDIST_ADJACENCIES_WEIGHT"
#define PIDIST_ADJACENCIES_WEIGHT_TYPE moab::MB_TYPE_INTEGER
#define PIDIST_PART_TAG_NAME "PIDIST_PART"
#define PIDIST_PART_TYPE moab::MB_TYPE_INTEGER
#define PIDIST_DISTRIBUTION_TIME_TAG_NAME "PIDIST_DISTRIBUTION_TIME"
#define PIDIST_DISTRIBUTION_TIME_TYPE moab::MB_TYPE_DOUBLE
#define PIDIST_COMMUNICATION_TIME_TAG_NAME "PIDIST_COMMUNICATION_TIME"
#define PIDIST_COMMUNICATION_TIME_TYPE moab::MB_TYPE_DOUBLE


typedef enum
{
	P_PartKway,
	P_PartGeomKway,
	P_PartGeom,
	P_RefineKway,
	P_AdaptiveRepart
} PIDist_ParmetisMethod;

typedef struct
{
	std::vector<double> wgt;
	int num;
} PIDist_wgt_pair;

typedef struct
{
	unsigned int from,to,id;
} PIDist_send_data;


int PIMisc_delete_shared_info(moab::Interface & imoab,EntityHandle set_handle, moab::ParallelComm & icomm)
{
	moab::Range ents;
	int result,dim;
	result = imoab.get_dimension(dim);
	if( result != MB_SUCCESS) RFAIL;
	for(int it = 0; it < dim; it++)
	{
		std::vector<unsigned char> pstat;
		std::vector<EntityHandle> sh;
		std::vector<int> sp;
		std::vector<EntityHandle> shs;
		std::vector<int> sps;
		ents.clear();
		result = imoab.get_entities_by_dimension(set_handle,it,ents);
		if( result != MB_SUCCESS) RFAIL;
		pstat.resize(ents.size(),0);
		sh.resize(ents.size(),0);
		sp.resize(ents.size(),-1);
		shs.resize(ents.size()*MAX_SHARING_PROCS,0);
		sps.resize(ents.size()*MAX_SHARING_PROCS,-1);
		result = imoab.tag_set_data(icomm.pstatus_tag(),ents,&pstat[0]);
		if( result != MB_SUCCESS) RFAIL;
		result = imoab.tag_set_data(icomm.sharedh_tag(),ents,&sh[0]);
		if( result != MB_SUCCESS) RFAIL;
		result = imoab.tag_set_data(icomm.sharedp_tag(),ents,&sp[0]);
		if( result != MB_SUCCESS) RFAIL;
		result = imoab.tag_set_data(icomm.sharedhs_tag(),ents,&shs[0]);
		if( result != MB_SUCCESS) RFAIL;
		result = imoab.tag_set_data(icomm.sharedps_tag(),ents,&sps[0]);
		if( result != MB_SUCCESS) RFAIL;
	}
	return MB_SUCCESS;
}


int PIMisc_print_part_data(moab::Interface & imoab,moab::ParallelComm & icomm, moab::Range ents, std::string comment)
{
	int result;
	unsigned int i;
	moab::Tag part_tag;
	std::vector <int> part_int(ents.size());
	result = imoab.tag_get_handle(PIDIST_PART_TAG_NAME,1,PIDIST_PART_TYPE,part_tag,moab::MB_TAG_DENSE|moab::MB_TAG_CREAT);
	if( result != MB_SUCCESS ) RFAIL;
	result = imoab.tag_get_data(part_tag,ents,&part_int[0]);
	if( result != MB_SUCCESS ) RFAIL;
	{
		std::map<int, int> part_num;
		std::map<int, int>::iterator part_num_i;
		for(i = 0; i < part_int.size(); i++)
		{
			part_num[part_int[i]]++;
		}
		{
			
			std::cout << comment << " on processor " << icomm.rank() << std::endl;
			for(part_num_i = part_num.begin(); part_num_i != part_num.end(); part_num_i++)
				std::cout << " have " << part_num_i->second << " elements of part " << part_num_i->first  << std::endl;
			std::cout << std::endl;
			
		}
	}
	return MB_SUCCESS;
}

int PIMisc_print_part_data(moab::Interface & imoab,moab::ParallelComm & icomm, std::vector<EntityHandle> ents, char * comment)
{
	int result;
	unsigned int i;
	moab::Tag part_tag;
	std::vector <int> part_int(ents.size());
	result = imoab.tag_get_handle(PIDIST_PART_TAG_NAME,1,PIDIST_PART_TYPE,part_tag,moab::MB_TAG_DENSE|moab::MB_TAG_CREAT);
	if( result != MB_SUCCESS ) RFAIL;
	result = imoab.tag_get_data(part_tag,&ents[0],ents.size(),&part_int[0]);
	if( result != MB_SUCCESS ) RFAIL;
	{
		std::map<int, int> part_num;
		std::map<int, int>::iterator part_num_i;
		for(i = 0; i < part_int.size(); i++)
		{
			part_num[part_int[i]]++;
		}
		{
			
			std::cout << "processor " << icomm.rank() << std::endl;
			std::cout << "comment " << comment << std::endl;
			for(part_num_i = part_num.begin(); part_num_i != part_num.end(); part_num_i++)
				std::cout << " have " << part_num_i->second << " elements of part " << part_num_i->first  << std::endl;
			
		}
	}
	return MB_SUCCESS;
}


int PIDist_distribute_entities(moab::Interface & imoab,EntityHandle set_handle, moab::ParallelComm & icomm, int num_local_parts = 1, PIDist_ParmetisMethod method = P_PartGeomKway)
{
	int debug_output = 1;
	int store_remote_handles = 0;
	int elem_dim = 3; //this parameter is used to look, which entities should be distributed
	int bridge_dim = 2; //this parameter is used to determine second-order adjacencies
	//these parameters are changed if global dimension is less then 3, look by ELEM_DIM_CHANGE
	
	
	//DONE//TASK: support for vertex weights, multiple per entity
	//DONE//TASK: support for adjacencies weights
	//DONE//TASK: automatically fill array of entities size
	//DONE//TASK: weight for every part set to support proccessors with different computing power
	//DONE//TASK: ITR should be selected automatically?
	// every call to exchange_tag should be measured by PIDist_add_communication_time()
	//DONE//TASK: should support multiple parts per proccess
	//DONE//TASK: replace std::find by red-black tree query
	//TASK: distribute part and weight information, while distributing graph information
	//DONE//TASK: reassemble "part" array in case AdaptiveRepart is requested and nparts != icomm.size()
	//TASK: tell to MOAB that we have part sets when MOAB will allow this

#if !defined(DEBUG_OUTPUT)
	debug_output = 0;
#endif

#if defined(STORE_REMOTE_HANDLES)
	store_remote_handles = 1;
#endif
	
	int result, dim, temp, have_vwgt = 0, have_adjwgt = 0;
	int part_div, part_offset;
	unsigned int i,j,k,nreserve = moab::MAX_SUB_ENTITIES, sent = 0,max_loc_parts;
	MPI_Comm comm = icomm.comm();
	idx_t nparts;
	//Array storing parts indexes
	std::vector<unsigned int> xparts;
	//These arrays determine local sparse matrix
	std::vector<idx_t> vtxdist;
	std::vector<idx_t> xadj;
	std::vector<idx_t> adjncy;
	//These arrays determine weights defined for vertices
	std::vector<idx_t> vwgt; //unused yet
	std::vector<idx_t> adjwgt; //unused yet
	std::vector<idx_t> vsize; //unused yet
	idx_t wgtflag = 0;
	idx_t numflag = 0;
	idx_t ndims;
	std::vector<real_t> xyz;
	//This sets weight for every part set in respect to according weight and tolerance
	idx_t ncon = 1;
	std::vector<real_t> tpwgts;
	std::vector<real_t> ubvec;
	//set options
	std::vector<idx_t> options(3,0);
	//number of edgecuts
	idx_t edgecut;
	//output data - which entity should belong to which part
	std::vector<idx_t> part;
	real_t itr = 1000.0, sum_tpwgts = 0;
	//send data to processors with zero nodes
	std::vector<PIDist_send_data> send;
	std::vector<int> adjncy_id(nreserve);
	idx_t mysize;
	double distribution_time;
	std::vector<double> local_tpwgtsd;
	std::vector<real_t> local_tpwgts;
	//to estimate memory used by entities
	std::vector<unsigned long> entity_storage;
	std::vector<int> recvcounts(icomm.size());
	std::vector<int> displs(icomm.size());
	std::vector<unsigned int> exchange_procs,part_data,procs_incoming_v;
	std::map<unsigned int,unsigned int>::iterator find_proc;
	moab::DataType tagtype;
	moab::Tag global_id, mesh_wgt_tag,ubvec_tag, vwgt_tag, adjwgt_tag,part_tag;
	moab::Tag dist_t_tag, comm_t_tag;
	moab::Range::iterator ent;
	moab::Range local,ghost;
	std::vector<moab::Range *> exchange_ents;
	std::set<unsigned int> part_set, procs_incoming, procs_outgoing;
	std::set<unsigned int>::iterator psiter;
	std::map<unsigned int,unsigned int> exchange_procs_helper,proc_by_part;
	std::vector<MPI_Request> sendReq, sendReqh;
	std::vector<MPI_Status> sendStat, sendStath;
	std::vector<moab::ParallelComm::Buffer> sendBuff, sendBuffh;
	std::vector<std::vector<EntityHandle> > L1hloc;
	std::vector<std::vector<EntityHandle> > L1hrem;
	std::vector<std::vector<int> > L1p;
	std::vector<EntityHandle> L2hloc;
	std::vector<EntityHandle> L2hrem;
	std::vector<unsigned int> L2p;				
	std::vector<EntityHandle> new_ents;
	moab::Range allsent,exch_tmp;
	moab::TupleList entprocs;
	moab::MeshTopoUtil topo(&imoab);
	
	
	distribution_time = MPI_Wtime();

	//compute mesh dimension
	result = imoab.get_dimension(dim);
	if( result != MB_SUCCESS ) RFAIL;
	
	//share dimension between processors
	temp = dim;
	result = MPI_Allreduce(&temp,&dim,1,MPI_INTEGER,MPI_MIN,icomm.comm());
	if( result != MPI_SUCCESS ) FAIL;
	
	if(!dim) dim = 3;
	ndims = dim;
	
	assert(dim > 0);
	
	if( dim < 3 )
	{
		//ELEM_DIM_CHANGE
		elem_dim = dim;
		bridge_dim = dim-1;
	}
	
	//compute part ranges for every processor
	//and total number of parts
	xparts.resize(icomm.size()+1);
	xparts[0] = 0;
	max_loc_parts = 0;
//	std::cout << icomm.rank() << " MPI_Allgather num_local_parts" <<std::endl;
	result = MPI_Allgather(&num_local_parts,1,MPI_UNSIGNED,(&xparts[0]+1),1,MPI_UNSIGNED,icomm.comm());
	if( result != MPI_SUCCESS ) FAIL;
	for(i = 1; i < xparts.size(); i++)
	{
		if( xparts[i] > max_loc_parts ) max_loc_parts = xparts[i];
		xparts[i] += xparts[i-1];
	}
	part_offset = (max_loc_parts*icomm.size()*icomm.rank());
	part_div = max_loc_parts*icomm.size();
		
	if( debug_output )
	{
		if( icomm.rank()==0 )
		{
			for(i = 0; i < xparts.size(); i++)
				std::cout << xparts[i] << " ";
			std::cout << std::endl;
		}
	}		

	nparts = xparts[xparts.size()-1];
	
	//compute processor # for every part
	for(i = 0; i < nparts; i++)
	{
		for(j = 0; j < icomm.size(); j++)
			if( i >= xparts[j] && i < xparts[j+1] )
			{
				proc_by_part[i] = j;
				break;
			}
	}
	
	//set graph size
	vtxdist.resize(icomm.size()+1);
	
	


	//set options
	if( debug_output )
	{	
		options[0] = 1;
		options[1] = 1;
		options[2] = 15;
	}
	
	
	
	//set ghost boundaries so we certanly know we can access adjacencies

	result = icomm.resolve_shared_ents(set_handle,elem_dim,bridge_dim);
	if( result != MB_SUCCESS ) RFAIL;
	result = icomm.exchange_ghost_cells(elem_dim,bridge_dim,1,0,1);
	if( result != MB_SUCCESS ) RFAIL;
	result = icomm.assign_global_ids(set_handle,elem_dim,0);
	if( result != MB_SUCCESS ) RFAIL;

	
	//find local and ghost entities
	result = PIMisc_split_entities(imoab,set_handle,icomm,local,ghost,elem_dim);
	if( result != MB_SUCCESS ) RFAIL;
	
	mysize = local.size();
	vtxdist[0] = 0;
//	std::cout << icomm.rank() << " MPI_Allgather mysize" <<std::endl;
	result = MPI_Allgather(&mysize,1,IDX_T,&vtxdist[1],1,IDX_T,icomm.comm());
	if( result != MPI_SUCCESS ) FAIL;
	for(i = 2; i < vtxdist.size(); i++)
		vtxdist[i] += vtxdist[i-1];


	if( !(method & P_PartGeom) )
	{	
		xadj.resize(local.size()+1);
		adjncy.reserve(local.size()*nreserve); //preallocate array
		xadj[0] = 0;
	}
	if( method & (P_PartGeom | P_PartGeomKway) )
		xyz.resize(local.size()*dim);
	part.resize(local.size());
				
	result = imoab.tag_get_handle(PARALLEL_GID_TAG_NAME,1,moab::MB_TYPE_INTEGER,global_id);
	if( result != MB_SUCCESS ) RFAIL;
	
	
	//compute itr parameter
	//if( method & P_AdaptiveRepart )
	{
		double comm_t, dist_t;
		j = 1;
		result = imoab.tag_get_handle(PIDIST_DISTRIBUTION_TIME_TAG_NAME,1,PIDIST_DISTRIBUTION_TIME_TYPE,dist_t_tag,moab::MB_TAG_SPARSE|moab::MB_TAG_CREAT);
		if( result == MB_TAG_NOT_FOUND ) j = 0;
		else if( result != MB_SUCCESS ) RFAIL;
		result = imoab.tag_get_handle(PIDIST_COMMUNICATION_TIME_TAG_NAME,1,PIDIST_COMMUNICATION_TIME_TYPE,comm_t_tag,moab::MB_TAG_SPARSE|moab::MB_TAG_CREAT);
		if( result == MB_TAG_NOT_FOUND ) j = 0;
		else if( result != MB_SUCCESS ) RFAIL;

		if( j )
		{
			result = imoab.tag_get_data(dist_t_tag,&set_handle,1,&dist_t);
			if( result == MB_TAG_NOT_FOUND ) j = 0;
			else if( result != MB_SUCCESS ) RFAIL;
			if( j )
			{
				result = imoab.tag_get_data(comm_t_tag,&set_handle,1,&comm_t);
				if( result == MB_TAG_NOT_FOUND ) j = 0;
				else if( result != MB_SUCCESS ) RFAIL;
				if( j ) itr = comm_t/dist_t;
			}
		}
		if( !j ) itr = 1000;
	}

	
	//estimate storage space
	if( method & P_AdaptiveRepart ) 
	{
		entity_storage.resize(local.size());
		imoab.estimated_memory_use(local,NULL,NULL,&entity_storage[0]);
		vsize.resize(local.size());
	}
	
	if( !(method & P_PartGeom ) )
	{
		result = imoab.tag_get_handle(PIDIST_VERTEX_WEIGHT_TAG_NAME,0,PIDIST_VERTEX_WEIGHT_TYPE,vwgt_tag,moab::MB_TAG_ANY);
		if( result == MB_SUCCESS ) have_vwgt = 1;
		else if( result == MB_TAG_NOT_FOUND ) have_vwgt = 0;
		else RFAIL;
		result = imoab.tag_get_handle(PIDIST_ADJACENCIES_WEIGHT_TAG_NAME,0,PIDIST_ADJACENCIES_WEIGHT_TYPE,adjwgt_tag,moab::MB_TAG_ANY);
		if( result == MB_SUCCESS ) have_adjwgt = 1;
		else if( result == MB_TAG_NOT_FOUND ) have_adjwgt = 0;
		else RFAIL;
		
		if( have_vwgt )
		{
			result = imoab.tag_get_data_type(vwgt_tag,tagtype);
			if( result != MB_SUCCESS ) RFAIL;
			if( tagtype == PIDIST_VERTEX_WEIGHT_TYPE )
			{
				result = imoab.tag_get_length(vwgt_tag,temp);
				if( result != MB_SUCCESS ) RFAIL;
				ncon = temp;
			}
			else have_vwgt = 0;
		}
		
		if( have_adjwgt )
		{
			result = imoab.tag_get_data_type(adjwgt_tag,tagtype);
			if( result != MB_SUCCESS ) RFAIL;
			if( tagtype == PIDIST_ADJACENCIES_WEIGHT_TYPE )
			{
				result = imoab.tag_get_length(adjwgt_tag,temp);
				if( result != MB_SUCCESS ) RFAIL;
				if( have_vwgt && ncon != temp ) FAIL;
			}
			else have_vwgt = 0;
		}
		if( have_vwgt ) vwgt.resize(local.size()*ncon);
		if( have_adjwgt) adjwgt.reserve(local.size()*ncon*nreserve);
		if( !have_adjwgt && !have_vwgt ) ncon = 1;
		//set constraint tolerances
		result = imoab.tag_get_handle(PIDIST_PART_WEIGHT_TOLERANCE_TAG_NAME,0,PIDIST_PART_WEIGHT_TOLERANCE_TYPE,ubvec_tag,moab::MB_TAG_MESH|moab::MB_TAG_ANY);
		j = 0;
		if( result == MB_SUCCESS )
		{
			result = imoab.tag_get_data_type(ubvec_tag,tagtype);
			if( result != MB_SUCCESS ) RFAIL;
			if( tagtype != PIDIST_PART_WEIGHT_TOLERANCE_TYPE) j = 1;
			else
			{
				result = imoab.tag_get_length(ubvec_tag,temp);
				if( result != MB_SUCCESS ) RFAIL;
				if( (have_adjwgt || have_vwgt) && temp != ncon ) j = 1;
				else
				{
					if( !have_adjwgt && !have_vwgt ) ncon = temp;
					ubvec.resize(ncon);
					result = imoab.tag_get_data(ubvec_tag,&set_handle,1,&ubvec[0]);
					if( result != MB_SUCCESS ) RFAIL;
				}
			}	
		}
		else if( result == MB_TAG_NOT_FOUND ) j = 1;
		else RFAIL;
		
		if(j) ubvec.resize(ncon,1.05);
		//set weights for every part
		local_tpwgts.resize(num_local_parts*ncon);
		result = imoab.tag_get_handle(PIDIST_PART_WEIGHT_TAG_NAME,0,PIDIST_PART_WEIGHT_TYPE,mesh_wgt_tag,moab::MB_TAG_MESH|moab::MB_TAG_ANY);
		j = 0;
		if( result == MB_SUCCESS )
		{
			result = imoab.tag_get_data_type(mesh_wgt_tag,tagtype);
			if( result != MB_SUCCESS ) RFAIL;
			if( tagtype != PIDIST_PART_WEIGHT_TYPE) j = 1;
			else
			{
				result = imoab.tag_get_length(mesh_wgt_tag,temp);
				if( result != MB_SUCCESS ) RFAIL;
				if( temp != num_local_parts*ncon ) j = 1;
				else
				{
					local_tpwgtsd.resize(ncon*num_local_parts);
					result = imoab.tag_get_data(mesh_wgt_tag,&set_handle,1,&local_tpwgtsd[0]);
					if( result != MB_SUCCESS ) RFAIL;
					for(i = 0; i < local_tpwgtsd.size(); i++) local_tpwgts[i] = local_tpwgtsd[i];
				}
			}
		}
		else if( result == MB_TAG_NOT_FOUND ) j = 1;
		else RFAIL;
		if( j )
		{
			for(i = 0; i < local_tpwgts.size(); i++)
				local_tpwgts[i] = 1.0/(double)nparts;
		}
		
		for(i = 0; i < icomm.size(); i++)
		{
			recvcounts[i] = (xparts[i+1]-xparts[i])*ncon;
			displs[i] = xparts[i]*ncon;
		}
		tpwgts.resize(nparts*ncon);
//		std::cout << icomm.rank() << " MPI_Allgather tpwgts" <<std::endl;
		result = MPI_Allgatherv(&local_tpwgts[0],num_local_parts*ncon,REAL_T,&tpwgts[0],&recvcounts[0],&displs[0],REAL_T,icomm.comm());
		if( result != MPI_SUCCESS ) FAIL;
		
		for(j = 0; j < ncon; j++)
		{
			sum_tpwgts = 0;
			for(i = 0; i < tpwgts.size()/ncon; i++) sum_tpwgts += tpwgts[i*ncon+j];
			for(i = 0; i < tpwgts.size()/ncon; i++) tpwgts[i*ncon+j] /= sum_tpwgts;
		}
	}
	
	//exchange weights
	if( have_vwgt )
	{
		result = icomm.exchange_tags(vwgt_tag,exch_tmp);
		if( result != MB_SUCCESS ) RFAIL;
	}
	if( have_adjwgt )
	{
		result = icomm.exchange_tags(adjwgt_tag,exch_tmp);
		if( result != MB_SUCCESS ) RFAIL;
	}
	
	//get part tag
	{
		result = imoab.tag_get_handle(PIDIST_PART_TAG_NAME,1,PIDIST_PART_TYPE,part_tag,moab::MB_TAG_DENSE|moab::MB_TAG_CREAT);
		if( result != MB_SUCCESS ) RFAIL;
	}
	
	if( method & P_AdaptiveRepart )
	{
		std::vector<int> part_int(local.size());
		result = imoab.tag_get_data(part_tag,local,&part_int[0]);
		for(i = 0; i < local.size(); i++)
			part[i] = part_int[i];
	}
	
	//make graph
	for(ent = local.begin(), i = 0; ent != local.end(); ent++, i++)
	{
		double position[3];
		moab::Range adjacencies_elements;
		if( method & (P_PartGeom | P_PartGeomKway) )
		{
			result = topo.get_average_position(*ent,position);
			if( result != MB_SUCCESS ) RFAIL;
			for(j = 0; j < (unsigned int)dim; j++)
				xyz[i*dim+j] = position[j];
		}
		if( !(method & P_PartGeom) )
		{
			if( !have_adjwgt )
			{
				result = topo.get_bridge_adjacencies(*ent,bridge_dim,elem_dim,adjacencies_elements);
				if( result != MB_SUCCESS ) RFAIL;
			}
			else
			{
				//here is the case when entity is shared by multiple edges
				//this happens with vertexes for example
				int r,q;
				std::map <EntityHandle,PIDist_wgt_pair> wgt_map;
				std::map <EntityHandle,PIDist_wgt_pair>::iterator iwgt_map;
				std::vector<double> oneadjwgt(ncon);
				moab::Range bridge_adjacencies;
				moab::Range::iterator bridge;
				result = imoab.get_adjacencies(&*ent,1,bridge_dim,0,bridge_adjacencies,moab::Interface::UNION);
				if( result != MB_SUCCESS ) RFAIL;
				for(bridge = bridge_adjacencies.begin(); bridge != bridge_adjacencies.end(); bridge++)
				{
					moab::Range adj;
					moab::Range::iterator iadj;
					
					result = imoab.tag_get_data(adjwgt_tag,&*bridge,1,&oneadjwgt[0]);
					if( result != MB_SUCCESS ) RFAIL;
					result = imoab.get_adjacencies(&*ent,1,elem_dim,0,adj,moab::Interface::UNION);
					if( result != MB_SUCCESS ) RFAIL;
					adj.erase(*ent);
					for(iadj = adj.begin(); iadj != adj.end(); iadj++)
					{					
						iwgt_map = wgt_map.find(*iadj);
						if( iwgt_map == wgt_map.end() )
						{
							wgt_map[*iadj].num = 1;
							iwgt_map = wgt_map.find(*iadj);
							iwgt_map->second.wgt.resize(ncon);
							for(q = 0; q < ncon; q++)
								iwgt_map->second.wgt[q] = oneadjwgt[q];
						}
						else
						{
							for(q = 0; q < ncon; q++)
								iwgt_map->second.wgt[q] += oneadjwgt[q];
							iwgt_map->second.num ++;
						}
					}
					adjacencies_elements.merge(adj);
				}
				for(q = 0, bridge = adjacencies_elements.begin(); bridge != adjacencies_elements.end(); q++, bridge++)
				{
					iwgt_map = wgt_map.find(*bridge);
					for(r = 0; r < ncon; r++)
						adjwgt.push_back(iwgt_map->second.wgt[r]/(double)iwgt_map->second.num);
				}
			}
			
			if( adjacencies_elements.size() > nreserve ) 
			{
				nreserve = adjacencies_elements.size();
				adjncy_id.reserve(nreserve); //reallocate array
			}
			
			adjncy_id.resize(adjacencies_elements.size());
			result = imoab.tag_get_data(global_id,adjacencies_elements,&adjncy_id[0]);
			if( result != MB_SUCCESS ) RFAIL;
			for(j = 0; j < adjacencies_elements.size(); j++) 
				adjncy.push_back(adjncy_id[j]);
			xadj[i+1] = xadj[i] + adjacencies_elements.size();
			
			//Fill weight arrays here
			if( have_vwgt )
			{
				result = imoab.tag_get_data(vwgt_tag,&*ent,1,&vwgt[i*ncon]);
				if( result != MB_SUCCESS ) RFAIL;
			}
		}
		if( method & P_AdaptiveRepart )
			vsize[i] = entity_storage[i];
	}
	
	
	
	wgtflag = have_vwgt + have_adjwgt*2;
	
	if( debug_output )
	{
		start_serial_output(icomm.comm());
		std::cout << "on processor " << icomm.rank() << std::endl;
		std::cout << "vtxdist: " << std::endl;
		for(i = 0; i < vtxdist.size(); i++)
		{
			std::cout << vtxdist[i] << " ";
			if( (i+1)%30 == 0 ) std::cout << std::endl;
		}
		std::cout << std::endl;
		if( !(method & P_PartGeom) )
		{
			std::cout << "xadj size " << xadj.size() << ":" << std::endl;
			for(i = 0; i < xadj.size(); i++)
			{
				std::cout << xadj[i] << " ";
				if( (i+1)%30 == 0 ) std::cout << std::endl;
			}
			std::cout << std::endl;
			std::cout << "adjncy size " << adjncy.size() << ":" << std::endl;
			for(i = 0; i < adjncy.size(); i++)
			{
				std::cout << adjncy[i] << " ";
				if( (i+1)%30 == 0 ) std::cout << std::endl;
			}
			std::cout << std::endl;
			if( have_vwgt )
			{
				std::cout << "vwgt: " << std::endl;
				for(i = 0; i < vwgt.size(); i++)
				{
					std::cout << vwgt[i] << " ";
					if( (i+1)%30 == 0 ) std::cout << std::endl;
				}
				std::cout << std::endl;
			}
			if( have_adjwgt )
			{
				std::cout << "adjwgt: " << std::endl;
				for(i = 0; i < adjwgt.size(); i++)
				{
					std::cout << adjwgt[i] << " ";
					if( (i+1)%30 == 0 ) std::cout << std::endl;
				}
				std::cout << std::endl;
			}
			if( method & P_AdaptiveRepart )
			{
				std::cout << "vsize: " << std::endl;
				for(i = 0; i < vsize.size(); i++)
				{
					std::cout << vsize[i] << " ";
					if( (i+1)%30 == 0 ) std::cout << std::endl;
				}
				std::cout << std::endl;
			}			
			std::cout << "wgtflag = " << wgtflag << std::endl;
			std::cout << "numflag = " << numflag << std::endl;
		}
		if( method & (P_PartGeom | P_PartGeomKway) )
		{
			std::cout << "ndims = " << ndims << std::endl;
			std::cout << "xyz: " << std::endl;
			for(i = 0; i < xyz.size(); i++)
			{
				std::cout << xyz[i] << " ";
				if( (i+1)%30 == 0 ) std::cout << std::endl;
			}
			std::cout << std::endl;
		}
		std::cout << "ncon = " << ncon << std::endl;
		std::cout << "nparts = " << nparts << std::endl;
		std::cout << "tpwgts: " << std::endl;
		for(i = 0; i < tpwgts.size(); i++)
		{
			std::cout << tpwgts[i] << " ";
			if( (i+1)%30 == 0 ) std::cout << std::endl;
		}
		std::cout << std::endl;
		std::cout << "ubvec: " << std::endl;
		for(i = 0; i < ubvec.size(); i++)
		{
			std::cout << ubvec[i] << " ";
			if( (i+1)%30 == 0 ) std::cout << std::endl;
		}
		std::cout << std::endl;
		std::cout << "options: " << std::endl;
		for(i = 0; i < options.size(); i++)
		{
			std::cout << options[i] << " ";
			if( (i+1)%30 == 0 ) std::cout << std::endl;
		}
		std::cout << std::endl;
		end_serial_output(icomm.comm());
	}
	
	//redistribute entities;
	{
		int it,jt,kt,flag;
		PIDist_send_data fill;
		send.reserve(nparts*2);
		for(it = icomm.size(); it >= 1; it--)
		{
			if( vtxdist[it]-vtxdist[it-1] == 0 )
			{
				fill.to = it - 1;
				flag = 0;
				for(jt = it-1; jt >= 1; jt--)
				{
					if( vtxdist[jt]-vtxdist[jt-1] > 0 )
					{
						fill.from = jt - 1;
						fill.id = vtxdist[jt]-1;
						send.push_back(fill);
						vtxdist[jt]--;
						for(kt = jt+1; kt < it; kt++)
							vtxdist[kt] = vtxdist[jt];
						flag = 1;
						break;
					}
				}
				if( !flag ) FAIL;
			}
		}
		
		for(i = 0; i < send.size(); i++)
		{
			idx_t trans_adjncy_size;
			real_t trans_xyz[3];
			std::vector<idx_t> trans_adjncy;
			if( send[i].to == send[i].from ) FAIL;
			if( send[i].from == icomm.rank() )
			{	
				if( !(method & P_PartGeom) )
				{
					trans_adjncy_size = xadj[send[i].id+1]-xadj[send[i].id];
					for(jt = xadj[send[i].id]; jt < xadj[send[i].id+1]; jt++)
						trans_adjncy.push_back(adjncy[jt]);
					result = MPI_Send(&trans_adjncy[0],trans_adjncy_size,IDX_T,send[i].to,2,icomm.comm());
					if( result != MPI_SUCCESS ) FAIL;
				}
				if( method & (P_PartGeom | P_PartGeomKway) )
				{
					for(jt = 0; jt < dim; jt++)
						trans_xyz[jt] = xyz[send[i].id*dim+jt];
					result = MPI_Send(trans_xyz,dim,REAL_T,send[i].to,3,icomm.comm());
					if( result != MPI_SUCCESS ) FAIL;
				}
			}
			else if( send[i].to == icomm.rank() )
			{
				MPI_Status stat;
				int msgsize;
				if( !(method & P_PartGeom) )
				{
					result = MPI_Probe(send[i].from,2,icomm.comm(),&stat);
					if( result != MPI_SUCCESS ) FAIL;
					result = MPI_Get_count(&stat,IDX_T,&msgsize);
					if( result != MPI_SUCCESS ) FAIL;
					trans_adjncy_size = msgsize;
					trans_adjncy.resize(trans_adjncy_size);
					result = MPI_Recv(&trans_adjncy[0],trans_adjncy_size,IDX_T,send[i].from,2,icomm.comm(),&stat);
					if( result != MPI_SUCCESS ) FAIL;
					xadj.resize(2);
					xadj[0] = 0;
					xadj[1] = trans_adjncy_size;
					for(jt = 0; jt < trans_adjncy_size; jt++)
						adjncy.push_back(trans_adjncy[jt]);
				}
				if( method & (P_PartGeom | P_PartGeomKway) )
				{
					result = MPI_Recv(&trans_xyz,dim,REAL_T,send[i].from,3,icomm.comm(),&stat);
					if( result != MPI_SUCCESS ) FAIL;
					xyz.resize(dim);
					for(jt = 0; jt < dim; jt++)
						xyz[jt] = trans_xyz[jt];
				}
				part.resize(1);
			}
		}
	}
	

	//select method and solve
	switch(method)
	{
		case P_PartKway:
		result = ParMETIS_V3_PartKway(&vtxdist[0],&xadj[0],&adjncy[0],&vwgt[0],&adjwgt[0],
									  &wgtflag,&numflag,&ncon,&nparts,&tpwgts[0],
									  &ubvec[0],&options[0],&edgecut,&part[0],&comm);
		break;
		case P_RefineKway:
		result = ParMETIS_V3_RefineKway(&vtxdist[0],&xadj[0],&adjncy[0],&vwgt[0],&adjwgt[0],
										&wgtflag,&numflag,&ncon,&nparts,&tpwgts[0],
										&ubvec[0],&options[0],&edgecut,&part[0],&comm);	
		break;
		case P_PartGeomKway:
		result = ParMETIS_V3_PartGeomKway(&vtxdist[0],&xadj[0],&adjncy[0],&vwgt[0],&adjwgt[0],
										  &wgtflag,&numflag,&ndims,&xyz[0],&ncon,&nparts,&tpwgts[0],
										  &ubvec[0],&options[0],&edgecut,&part[0],&comm);
		break;
		case P_PartGeom:
		result = ParMETIS_V3_PartGeom(&vtxdist[0],&ndims,&xyz[0],&part[0],&comm);
		break;
		case P_AdaptiveRepart:
		result = ParMETIS_V3_AdaptiveRepart(&vtxdist[0],&xadj[0],&adjncy[0],&vwgt[0],&vsize[0],&adjwgt[0],
											&wgtflag,&numflag,&ncon,&nparts,&tpwgts[0],
											&ubvec[0],&itr,&options[0],&edgecut,&part[0],&comm);
		break;
	}
	
	if( result != METIS_OK ) FAIL;
	
	//distribute computed parts back
	{
		for(i = 0; i < send.size(); i++)
		{
			if( send[i].to == icomm.rank() )
			{
				result = MPI_Send(&part[0],1,IDX_T,send[i].from,4,icomm.comm());
				if( result != MPI_SUCCESS ) FAIL;
				part.clear();
			}
			else if( send[i].from == icomm.rank() )
			{ 
				result = MPI_Recv(&part[send[i].id],1,IDX_T,send[i].to,4,icomm.comm(),MPI_STATUS_IGNORE);
				if( result != MPI_SUCCESS ) FAIL;
			}
		}
	}
							 
	//debug
	if( debug_output )
	{
		start_serial_output(icomm.comm());
		std::cout << "distribution on " << icomm.rank() << std::endl;
		for(i = 0; i < part.size(); i++)
		{
			std::cout << part[i] << " ";
			if( (i+1)%30 == 0 ) std::cout << std::endl;
		}
		std::cout << std::endl;
		end_serial_output(icomm.comm());
	}
	
	//assign entities to part sets
	{
		std::vector<int> part_int(part.size());
		for(i = 0; i < part.size(); i++) part_int[i] = part[i];
		result = imoab.tag_set_data(part_tag,local,&part_int[0]);
		if( result != MB_SUCCESS ) RFAIL;
		result = icomm.exchange_tags(part_tag,exch_tmp);
		if( result != MB_SUCCESS ) RFAIL;
		
		if( debug_output )
		{
			start_serial_output(icomm.comm());
			PIMisc_print_part_data(imoab,icomm,local,"local without ghost");
			end_serial_output(icomm.comm());
		}
	}
	
	if( debug_output )
	{
		moab::Range ents;
		result = imoab.get_entities_by_dimension(set_handle,elem_dim,ents);
		if( result != MB_SUCCESS ) RFAIL;
		start_serial_output(icomm.comm());
		PIMisc_print_part_data(imoab,icomm,ents,"local with ghost");
		end_serial_output(icomm.comm());
	}
	
	
	//compute exchange arrays of processors and entities
	{
		for(i = 0; i < part.size(); i++)
		{
			if( proc_by_part[part[i]] == icomm.rank() ) continue;
			j = proc_by_part[part[i]];
			find_proc = exchange_procs_helper.find(j);
			if( find_proc == exchange_procs_helper.end() )
			{
				k = exchange_procs.size();
				exchange_procs.push_back(j);
				exchange_ents.push_back(new moab::Range);
				exchange_procs_helper[j] = k;
			} else k = find_proc->second;
			exchange_ents[k]->insert(local[i]);
		}
		
		if( debug_output )
		{
			start_serial_output(icomm.comm());
			std::cout << " on processor " << icomm.rank() << std::endl;
			for(i = 0; i < exchange_procs.size(); i++)
			{
				std::cout << " send " << exchange_ents[i]->subset_by_dimension(3).size()
						  << " entities and " << exchange_ents[i]->subset_by_dimension(0).size()
						  << " vertexes to processor " << exchange_procs[i] << std::endl;
			}
			end_serial_output(icomm.comm());
		}
		
		for(i = 0; i < exchange_ents.size(); i++)
		{
			std::pair<Range::const_iterator, Range::const_iterator>
			set_range = exchange_ents[i]->equal_range(MBENTITYSET);
			for (Range::const_iterator rit = set_range.first; rit != set_range.second; rit++) 
			{
				result = imoab.get_entities_by_type(*rit, MBVERTEX, *exchange_ents[i]);
				if( result != MB_SUCCESS ) RFAIL;
			}
			if( result != MB_SUCCESS ) RFAIL;
			Range tmp_ents;
			std::copy(exchange_ents[i]->begin(), set_range.first, range_inserter(tmp_ents));
			result = imoab.get_adjacencies(tmp_ents, 0, 0, *exchange_ents[i],moab::Interface::UNION);
			if( result != MB_SUCCESS ) RFAIL;
			
			if( store_remote_handles )
			{
				moab::Range tmp_range;
				result = icomm.filter_pstatus(*exchange_ents[i], PSTATUS_SHARED, PSTATUS_AND,exchange_procs[i], &tmp_range);
				if( result != MB_SUCCESS ) RFAIL;
				if( !tmp_range.empty() )
					*exchange_ents[i] = subtract(*exchange_ents[i],tmp_range);
				if( exchange_ents[i]->empty() )
				{
					delete exchange_ents[i];
					exchange_ents.erase(exchange_ents.begin()+i);
					exchange_procs.erase(exchange_procs.begin()+i);
					i--;
				}
			}
		}
	}
	
	
	
	if( debug_output )
	{
		start_serial_output(icomm.comm());
		std::cout << " on processor " << icomm.rank() << std::endl;
		for(i = 0; i < exchange_procs.size(); i++)
		{
			std::cout << " send " << exchange_ents[i]->subset_by_dimension(3).size()
					  << " entities and " << exchange_ents[i]->subset_by_dimension(0).size()
					  << " vertexes to processor " << exchange_procs[i] << std::endl;
		}
		end_serial_output(icomm.comm());
	}

#if defined(USE_EXCHANGE_OWNED_MESHES) || defined(USE_SENDRECV)
	//compute if there are processors that want to communicate with us	
	{
		std::vector<int> incoming(icomm.size(),0);
		std::vector<int> incoming_mat(icomm.size()*icomm.size(),0);
		for(i = 0; i <exchange_procs.size(); i++)
			incoming[exchange_procs[i]] = 1;
		std::cout << icomm.rank() << " MPI_Allgather incoming" <<std::endl;
		MPI_Allgather(&incoming[0],icomm.size(),MPI_INT,&incoming_mat[0],icomm.size(),MPI_INT,icomm.comm());
		for(i = 0; i < exchange_procs.size(); i++)
			procs_outgoing.insert(exchange_procs[i]);
		for(i = 0; i < incoming.size(); i++)
			if( i != icomm.rank() && incoming_mat[i*icomm.size()+icomm.rank()] == 1 )
				procs_incoming.insert(i);
		for(psiter = procs_incoming.begin(); psiter != procs_incoming.end(); psiter++)
			procs_incoming_v.push_back(*psiter);
	}
#endif
	
	//now send entities to other processors
#if defined(USE_EXCHANGE_OWNED_MESHES)
	{
		std::vector<MPI_Request> recv_ent_reqs, recv_remoteh_reqs;
//		if( procs_incoming_v.size() )
//		{
//			result = icomm.post_irecv(procs_incoming_v);
//			if( result != MB_SUCCESS ) RFAIL;
//		}
		result = icomm.exchange_owned_meshs(exchange_procs,exchange_ents,recv_ent_reqs,recv_remoteh_reqs,1,0,1,3);
		if( result != MB_SUCCESS ) RFAIL;
	}
#elif defined(USE_SENDRECV)
	{
		int incoming1 = 0, incoming2 = 0;
		if( exchange_procs.size() )
		{
			result = icomm.send_entities(exchange_procs,exchange_ents,incoming1,incoming2,1);
			if( result != MB_SUCCESS ) RFAIL;
		}
		//std::cout << icomm.rank() << " incoming1: " << incoming1 << " incoming2: " << incoming2 << std::endl;
		{
			//post irecv for incoming remote handles??
			if( procs_incoming_v.size() )
			{
				result = icomm.post_irecv(exchange_procs);
				if( result != MB_SUCCESS ) RFAIL;
			}
			result = icomm.recv_entities(procs_incoming,procs_incoming_v.size()+incoming1,incoming2,1,1);
			if( result != MB_SUCCESS ) RFAIL;
		}
	}
#else
	//Following code rely on inner functionality of moab
	{
		//fill tuples
		{
			unsigned int npairs = 0;
			for (i = 0; i < exchange_procs.size(); i++) 
			{
				int n_ents = exchange_ents[i]->size();
				if (n_ents > 0) 
				{
					npairs += n_ents; // get the total # of proc/handle pairs
					allsent.merge(*exchange_ents[i]);
				}
			}
			entprocs.initialize(1, 0, 1, 0, npairs);
			entprocs.enableWriteAccess();
			// put the proc/handle pairs in the list
			for (i = 0; i < exchange_procs.size(); i++) 
			{
				for (moab::Range::iterator rit = exchange_ents[i]->begin(); rit != exchange_ents[i]->end(); rit++) 
				{
					entprocs.vi_wr[entprocs.get_n()] = exchange_procs[i];
					entprocs.vul_wr[entprocs.get_n()] = *rit;
					entprocs.inc_n();
				}
			}		
			// sort by handle
			moab::TupleList::buffer sort_buffer;
			sort_buffer.buffer_init(npairs);
			entprocs.sort(1, &sort_buffer);
			entprocs.disableWriteAccess();
			sort_buffer.reset();
			allsent.compactness();
		}
		
		
		//compute global pairs of send_to/recv_from
		{
			PIDist_send_data fill;
			std::vector<unsigned int> exchange_procs_global;
			unsigned int max_loc_size, loc_size = exchange_procs.size(), nreq = 0, nreqh = 0;
			send.resize(0);
			result = MPI_Allreduce(&loc_size,&max_loc_size,1,MPI_UNSIGNED,MPI_MAX,icomm.comm());
			if( result != MPI_SUCCESS ) FAIL;
			exchange_procs.resize(max_loc_size,UINT_MAX);
			exchange_procs_global.resize(max_loc_size*icomm.size(),UINT_MAX);
//			std::cout << icomm.rank() << " MPI_Allgather exchange_procs" <<std::endl;
			result = MPI_Allgather(&exchange_procs[0],max_loc_size,MPI_UNSIGNED,&exchange_procs_global[0],max_loc_size,MPI_UNSIGNED,icomm.comm());
			if( result != MPI_SUCCESS ) FAIL;
			for(i = 0; i < icomm.size(); i++)
				for(j = 0; (j < max_loc_size) && (exchange_procs_global[i*max_loc_size+j] != UINT_MAX); j++)
				{
					fill.from = i;
					fill.to = exchange_procs_global[i*max_loc_size+j];
					if( fill.from == icomm.rank() )
						fill.id = j;
					send.push_back(fill);
					if( fill.from == icomm.rank() ) nreq++;
					if( fill.to == icomm.rank() ) nreqh++;
				}
			sendReq.resize(nreq,MPI_REQUEST_NULL);
			sendStat.resize(nreq);
			sendBuff.resize(nreq,moab::ParallelComm::Buffer(1000000));
			
			L1hloc.resize(nreqh);
			L1hrem.resize(nreqh);
			L1p.resize(nreqh);
			
			L2hloc.resize(nreqh);
			L2hrem.resize(nreqh);
			L2p.resize(nreqh);
			
			sendReqh.resize(nreqh,MPI_REQUEST_NULL);
			sendStath.resize(nreqh);
			sendBuffh.resize(nreqh,moab::ParallelComm::Buffer(1000000));

		}
		
		if( debug_output )
		{
			start_serial_output(icomm.comm());
			std::cout << icomm.rank() <<  " Started data exchange" << std::endl;
			end_serial_output(icomm.comm());
		}
		{

			sent = 0;
			for(i = 0; i < send.size(); i++)
			{
				
				if( send[i].to == send[i].from ) continue; //don't touch already owned entities
				if( send[i].from == icomm.rank() ) //send buffer
				{
					
					int buff_size;
					//int ind = 
					//icomm.get_buffers(send[i].to); //let moab know about processor we communicate with
					sendBuff[sent].reset_buffer(sizeof(int));
					result = icomm.pack_buffer(*exchange_ents[send[i].id],
											   0,//adj
											   1,//tags
											   store_remote_handles,//remoteh
											   send[i].to,&(sendBuff[sent]),&entprocs,&allsent);
					if( result != MB_SUCCESS ) RFAIL;
					sendBuff[sent].set_stored_size();
					buff_size = sendBuff[sent].get_stored_size();
					if( debug_output )
					{
						std::cout << " MPI_Isend " << send[i].from << "->" 
								  << send[i].to << " size "
								  << buff_size << " id " << send[i].id << std::endl;
					}
					result = MPI_Isend(sendBuff[sent].mem_ptr,buff_size,MPI_UNSIGNED_CHAR,
									   send[i].to,i,icomm.comm(),&sendReq[sent]);
					if( result != MPI_SUCCESS ) FAIL;
					sent++;
					
					
					
				}
			}
			//imoab.remove_entities(icomm.partition_sets().front(),allsent);
			//if( result != MB_SUCCESS ) RFAIL;
			entprocs.reset();
			for(i = 0; i < send.size(); i++)
			{
				if( send[i].to == send[i].from ) continue; //don't touch already owned entities
				if( send[i].to == icomm.rank() ) //recieve buffer
				{
					std::vector<EntityHandle> new_ents;
					moab::ParallelComm::Buffer buff(1000000);
					MPI_Status stat;
					//int ind = icomm.get_buffers(send[i].from);
					int buff_size, tmp;
					result = MPI_Probe(send[i].from,i,icomm.comm(),&stat);
					if( result != MPI_SUCCESS ) FAIL;
					result = MPI_Get_count(&stat,MPI_UNSIGNED_CHAR,&buff_size);
					if( result != MPI_SUCCESS ) FAIL;
					buff.reserve(buff_size);
					if( debug_output )
					{
						std::cout << " MPI_Recv " << send[i].to << "<-" 
								  << send[i].from << " size " 
								  << buff_size << std::endl;
					}
					result = MPI_Recv(buff.mem_ptr,buff_size,MPI_UNSIGNED_CHAR,
									  send[i].from,i,icomm.comm(),&stat);
					if( result != MPI_SUCCESS ) FAIL;				
					result = MPI_Get_count(&stat,MPI_UNSIGNED_CHAR,&tmp);
					if( result != MPI_SUCCESS ) FAIL;
					if( tmp != buff_size ) FAIL;
					buff.reset_ptr(sizeof(int));
					result = icomm.unpack_buffer(buff.buff_ptr,
												 store_remote_handles,//remoteh
												 send[i].from,
												 -1,//ind,//ind
												 L1hloc,L1hrem,L1p,
												 L2hloc,L2hrem,L2p,new_ents,
												 1//created_iface
												 );
					if( result != MB_SUCCESS ) RFAIL;
					if( debug_output )
						std::cout << " proc " << icomm.rank() << " unpacked " << new_ents.size() << " entities " << std::endl;
					//imoab.add_entities(icomm.partition_sets().front(), &new_ents[0], new_ents.size());
					//if( result != MB_SUCCESS ) RFAIL;
				}
			}
			//imoab.add_entities(icomm.partition_sets().front(), &new_ents[0], new_ents.size());
			//if( result != MB_SUCCESS ) RFAIL;
			
			
			//bother about messages we have sent
			result = MPI_Waitall(sendReq.size(),&sendReq[0],&sendStat[0]);
			if( result != MPI_SUCCESS ) FAIL;
			
			if( store_remote_handles ) //exchange remote handles
			{
				sent = 0;
				for(i = 0; i < send.size(); i++)
				{
					if( send[i].to == send[i].from ) continue; //don't touch already owned entities
					if( send[i].to == icomm.rank() ) // send remote handles
					{
						int buff_size;
						int ind = icomm.get_buffers(send[i].from);
						sendBuffh[sent].reset_buffer(sizeof(int));
						//if( ind >= L1hloc.size() ) std::cout << "ERROR HERE" << std::endl;
						result = icomm.pack_remote_handles(L1hloc[ind],L1hrem[ind],L1p[ind],send[i].from,&sendBuffh[sent]);
						if( result != MB_SUCCESS ) RFAIL;
						sendBuffh[sent].set_stored_size();
						buff_size = sendBuffh[sent].get_stored_size();
					
						if( debug_output )
						{
							std::cout << " MPI_Isend handles " << send[i].to << "->" 
									  << send[i].from << " size "
									  << buff_size << std::endl;
						}
					
						result = MPI_Isend(sendBuffh[sent].mem_ptr,buff_size,MPI_UNSIGNED_CHAR,
										   send[i].from,send.size()+i,icomm.comm(),&sendReqh[sent]);
						if( result != MPI_SUCCESS ) FAIL;
						sent++;
						//exchange remote handles
					}
				}
				for(i = 0; i < send.size(); i++)
				{
					if( send[i].to == send[i].from ) continue; //don't touch already owned entities
					if( send[i].from == icomm.rank() ) //recv remote handles
					{
						moab::ParallelComm::Buffer buff(1000000);
						MPI_Status stat;
						int buff_size, tmp;
						result = MPI_Probe(send[i].to,i+send.size(),icomm.comm(),&stat);
						if( result != MPI_SUCCESS ) FAIL;
						result = MPI_Get_count(&stat,MPI_UNSIGNED_CHAR,&buff_size);
						if( result != MPI_SUCCESS ) FAIL;
						buff.reserve(buff_size);
						
						if( debug_output )
						{
							std::cout << " MPI_Recv handles " << send[i].from << "<-" 
									  << send[i].to << " size " 
									  << buff_size << std::endl;
						}
						
						result = MPI_Recv(buff.mem_ptr,buff_size,MPI_UNSIGNED_CHAR,
										  send[i].to,i+send.size(),icomm.comm(),&stat);
						if( result != MPI_SUCCESS ) FAIL;	
						result = MPI_Get_count(&stat,MPI_UNSIGNED_CHAR,&tmp);
						if( result != MPI_SUCCESS ) FAIL;
						if( tmp != buff_size ) FAIL;
						buff.reset_ptr(sizeof(int));
						result = icomm.unpack_remote_handles(send[i].to,buff.buff_ptr,L2hloc,L2hrem,L2p);
						if( result != MB_SUCCESS ) RFAIL;
						
					}
					
					//MPI_Barrier(icomm.comm());
				}
				result = MPI_Waitall(sendReqh.size(),&sendReqh[0],&sendStath[0]);
				if( result != MPI_SUCCESS ) FAIL;
			}
			
			if( debug_output )
			{
				std::cout << icomm.rank() << " DONE UNPACKING " << std::endl;
			}
		}
	}
#endif	

	//delete allocated ranges
	{
		for(i = 0; i < exchange_ents.size(); i++)
			delete exchange_ents[i];
	}
	
	
	//look for parts, delete parts that we shouldn't own
	//merge duplicate parts
	//start_serial_output(icomm.comm());
	
	{
		std::vector<int> part_int;
		std::map<int, int> part_num;
		std::map<int, int>::iterator part_num_i;
		local.clear();
		int erase = 0;
		result = imoab.get_entities_by_dimension(set_handle,elem_dim,local);
		if( result != MB_SUCCESS ) RFAIL;
		part_int.resize(local.size());
		result = imoab.tag_get_data(part_tag,local,&part_int[0]);
		if( result != MB_SUCCESS ) RFAIL;
		for(i = 0; i < part_int.size(); i++)
		{
			if( !((unsigned int)part_int[i] >= xparts[icomm.rank()] && (unsigned int)part_int[i] < xparts[icomm.rank()+1]) )
			{
				EntityHandle del = local[i];
				result = imoab.delete_entities(&del,1);
				if( result != MB_SUCCESS ) RFAIL;
				erase++;
			}
			part_num[part_int[i]]++;
		}
		if( debug_output )
		{
			start_serial_output(icomm.comm());
			std::cout << icomm.rank() << " deleted " << erase << " of " << local.size() << " entities " << std::endl;
			for(part_num_i = part_num.begin(); part_num_i != part_num.end(); part_num_i++)
				std::cout << " have " << part_num_i->second << " elements of part " << part_num_i->first << std::endl;
			end_serial_output(icomm.comm());
		}

	}
	
	//end_serial_output(icomm.comm());
	
	if( debug_output )
		std::cout << icomm.rank() << " finished with sets" << std::endl;
	
	//now we store only local entities, delete non-local adjacencies
	//look for adjacencies we have to delete
	
	{
		moab::Range sets;
		local.clear();
		result = imoab.get_entities_by_dimension(set_handle,elem_dim,local);
		if( result != MB_SUCCESS ) RFAIL;
		result = imoab.get_entities_by_type(set_handle,moab::MBENTITYSET,sets);
		if( result != MB_SUCCESS ) RFAIL;
		for(int it = dim; it >= 0; it--)
		{
			if(it != elem_dim )
			{
				moab::Range erase, keep, sets;
				result = imoab.get_entities_by_dimension(set_handle,it,erase);
				if( result != MB_SUCCESS ) RFAIL;
				result = imoab.get_adjacencies(local,it,1,keep,moab::Interface::UNION);
				if( result != MB_SUCCESS ) RFAIL;				
				if( !keep.empty() )
					erase -= keep;
				result = imoab.delete_entities(erase);
				if( result != MB_SUCCESS ) RFAIL;
				
				if( debug_output )
					std::cout << icomm.rank() << " erasing adjacencies " << it 
							  << " size " << erase.size() << std::endl;
			}
		}
		if(!store_remote_handles )PIMisc_delete_ghost(imoab,set_handle,icomm,elem_dim);
	}
	
	
	if( debug_output )
		std::cout << icomm.rank() << " finished with adjacencies" << std::endl;
	
	
	//now compute partition set
	{
		unsigned int rank = icomm.rank();
		moab::Tag parallel_partition = icomm.part_tag();
		moab::Range parallel_sets;
		result = imoab.get_entities_by_type_and_tag(set_handle,moab::MBENTITYSET,&parallel_partition,NULL,1,parallel_sets,moab::Interface::UNION);
		if( result != MB_SUCCESS ) RFAIL;
		if( !parallel_sets.empty() )
		{
			result = imoab.clear_meshset(parallel_sets);
			if( result != MB_SUCCESS ) RFAIL;
		}
		if( parallel_sets.size() == 0 )
		{
			EntityHandle new_set;
			result = imoab.create_meshset(moab::MESHSET_SET,new_set);
			if( result != MB_SUCCESS ) RFAIL;
			parallel_sets.insert(new_set);
		}
		else if( parallel_sets.size() != 1 )
		{
			unsigned int num_del = parallel_sets.size()-1;
			for(i = 0; i < num_del; i++)
			{
				EntityHandle old_set = parallel_sets.pop_back();
				result = imoab.delete_entities(&old_set,1);
				if( result != MB_SUCCESS ) RFAIL;
			}
		}
		result = imoab.tag_set_data(parallel_partition,parallel_sets,&rank);
		if( result != MB_SUCCESS ) RFAIL;
		local.clear();
		result = imoab.get_entities_by_dimension(set_handle,elem_dim,local);
		if( result != MB_SUCCESS ) RFAIL;
		result = imoab.add_entities(parallel_sets[0],local);
		if( result != MB_SUCCESS ) RFAIL;
		result = imoab.get_entities_by_type_and_tag(set_handle,MBENTITYSET,&parallel_partition,NULL,1,icomm.partition_sets());
		if( result != MB_SUCCESS ) RFAIL;
	}
	
	//fill information about shared ents
	{
		result = PIMisc_delete_shared_info(imoab,set_handle,icomm);
		//if( result != MB_SUCCESS ) RFAIL;
		//result = icomm.create_interface_sets(elem_dim,bridge_dim);
		//if( result != MB_SUCCESS ) RFAIL;
		icomm.sharedEnts.clear();
		icomm.interface_sets().clear();
		icomm.reset_all_buffers();
		/*
		std::map<std::vector<int>, std::vector<EntityHandle> > proc_nvecs;
		int procs[MAX_SHARING_PROCS];
		EntityHandle handles[MAX_SHARING_PROCS];
		int nprocs;
		unsigned char pstat;
		for (std::vector<EntityHandle>::iterator vit = icomm.sharedEnts.begin(); vit != icomm.sharedEnts.end(); vit++) 
		{
			if (imoab.dimension_from_handle(*vit) > 2)
				continue;
			result = icomm.get_sharing_data(*vit, procs, handles, pstat, nprocs);
			if( result != MB_SUCCESS ) RFAIL;
			std::sort(procs, procs+nprocs);
			std::vector<int> tmp_procs(procs, procs + nprocs);
			assert(tmp_procs.size() != 2);
			proc_nvecs[tmp_procs].push_back(*vit);
		}
		
		// create interface sets from shared entities
		result = icomm.create_interface_sets(proc_nvecs, elem_dim, bridge_dim);
		if( result != MB_SUCCESS ) RFAIL;
		*/
		//PIMisc_delete_shared_info(imoab,set_handle,icomm);
		
	}
	
	
	if( debug_output )
	{
		//start_serial_output(icomm.comm());
		std::cout << icomm.rank() << " distribution done" << std::endl;
		//end_serial_output(icomm.comm());
	}
	
	distribution_time = MPI_Wtime() - distribution_time;
	
	//attach distribution time to tag, reset communication time
	{
		double dist_t,comm_t = 0;
		result = MPI_Allreduce(&distribution_time,&dist_t,1,MPI_DOUBLE,MPI_MAX,icomm.comm());
		if( result != MPI_SUCCESS ) FAIL;
		result = imoab.tag_set_data(dist_t_tag,&set_handle,1,&dist_t);
		if( result != MB_SUCCESS ) RFAIL;
		result = imoab.tag_set_data(comm_t_tag,&set_handle,1,&comm_t);
		if( result != MB_SUCCESS ) RFAIL;

	}
	
	return MB_SUCCESS;
}

int print_ents_by_dim(moab::Interface & imoab,const EntityHandle set_handle, moab::ParallelComm & icomm)
{
	int num,dims,result;
	start_serial_output(icomm.comm());
	result = imoab.get_dimension(dims);
	if( result != MB_SUCCESS ) {end_serial_output(icomm.comm()); RFAIL;}
	std::cout << "on processor " << icomm.rank() << " mesh dimension " << dims << std::endl;
	for(int dim = 0; dim <= dims; dim++)
	{
		result = imoab.get_number_entities_by_dimension(set_handle,dim,num);
		if( result != MB_SUCCESS ) {end_serial_output(icomm.comm()); RFAIL;}
		std::cout << "elements of dimension " << dim << ": " << num << std::endl;
	}
	if( icomm.rank() != icomm.size()-1 ) std::cout << std::endl;
	end_serial_output(icomm.comm());
	return MB_SUCCESS;
}


double dist(double v1[3], double v2[3])
{
	return sqrt((v2[0]-v1[0])*(v2[0]-v1[0])+(v2[1]-v1[1])*(v2[1]-v1[1])+(v2[2]-v1[2])*(v2[2]-v1[2]));
}

int dimdir(double v1[3], double v2[3])
{
	for(int i = 0; i < 3; i++)
		if( fabs(v1[i]-v2[i]) > 1e-8 )
			return i;
	return -1;
}

//Starting point
int main(int argc,char ** argv)
{
	int rank,size,result;
	
	MPI_Init(&argc,&argv);
	MPI_Comm_rank(MPI_COMM_WORLD,&rank);
	MPI_Comm_size(MPI_COMM_WORLD,&size);
	
	if( argc > 1 )
	{
		moab::Range exch_tmp;
		moab::Core moab;
		moab::Interface & imoab = moab;
		moab::Tag phi_tag, global_id;
		moab::ParallelComm pcomm(&imoab,MPI_COMM_WORLD);
		moab::Range local,ghost;
		std::string version;
		
		imoab.impl_version(&version);
		
		
		if( pcomm.rank() == 0 ) std::cout << "MOAB version: " << version << std::endl;
		
		pcomm.set_debug_verbosity(0);
		
		//Loading mesh into program
		{
			if( argc > 2 )
			{
				switch(atoi(argv[2]))
				{
					case 1: result = PIFile_load_distributed(imoab,imoab.get_root_set(),pcomm,argv[1]); break;
					case 2: result = PIFile_load_serial(imoab,imoab.get_root_set(),pcomm,argv[1]); break;
					default:
					case 0: result = PIFile_load_serial_rank(imoab,imoab.get_root_set(),pcomm,argv[1],0); break;
				}
			}
			else result = PIFile_load_serial_rank(imoab,imoab.get_root_set(),pcomm,argv[1],0);
			CHKERR("load failed on files %s\n",argv[1]);
		}
		
		
		print_ents_by_dim(imoab,imoab.get_root_set(),pcomm);
		
		result = PIDist_distribute_entities(imoab,imoab.get_root_set(),pcomm,1,P_PartGeomKway);
		CHKERR("cannot distribute ents\n");
		
		
		print_ents_by_dim(imoab,imoab.get_root_set(),pcomm);
		
		{
			for(int i = 0; i <= 3; i++)
			{
				local.clear(), ghost.clear();
				PIMisc_split_entities(imoab,imoab.get_root_set(),pcomm,local,ghost,i);
				start_serial_output(pcomm.comm());
				std::cout << "rank " << rank << " have " << local.size() 
						  << " local elements and " << ghost.size() << " ghost elements " 
						  << " of dimension " << i << std::endl;
				end_serial_output(pcomm.comm());
			}
		}
		
		
		
		//prepare ghost boundaries
		{
			result = pcomm.resolve_shared_ents(imoab.get_root_set(),3,2);
			CHKERR("cannot resolve ents\n");
			
			result = pcomm.exchange_ghost_cells(3,2,1,0,1);
			CHKERR("cannot exchange ghost\n");
			
			//return 0;
		}

		
		
		//balance mesh
		/*
		{
			MBZoltan zoltan(&imoab,false,argc,argv);
			result = zoltan.balance_mesh("parmetis","refinekway");
		}
		*/
		
		//gather array of local and ghost entities for traversal
		{
			local.clear(),ghost.clear();
			PIMisc_split_entities(imoab,imoab.get_root_set(),pcomm,local,ghost,3);
			start_serial_output(pcomm.comm());
			std::cout << "rank " << rank << " have " << local.size() << " local elements and " << ghost.size() << " ghost elements " << std::endl;
			end_serial_output(pcomm.comm());	
		}
		//return 0;
		
		/*
		//Attach tag to entities
		{
			
			std::vector<double> values(local.size(),0);
			result = imoab.tag_get_handle("TAG_PHI",1,moab::MB_TYPE_DOUBLE,phi_tag,moab::MB_TAG_DENSE|moab::MB_TAG_CREAT);
			CHKERR("cannot create tag\n");
			result = imoab.tag_set_data(phi_tag,local,&values[0]);
			CHKERR("cannot set tag data to entities\n");
			result = pcomm.exchange_tags(phi_tag,exch_tmp);
			CHKERR("cannot exchange ghost data\n");
			result = pcomm.assign_global_ids(imoab.get_root_set(),3,0);
			CHKERR("cannot assign global id\n");
			result = imoab.tag_get_handle(PARALLEL_GID_TAG_NAME,1,moab::MB_TYPE_INTEGER,global_id);
			CHKERR("cannot get global id tag\n");
		}
		
		
		//Assemble matrix in petsc
		//solve \nabla \cdot \nabla phi = 0 equation
		//with boundary conditions
		//and put solution back to mesh
		
		{
			int its; double norm;
			KSP ksp = PETSC_NULL;
			Vec x = PETSC_NULL,b = PETSC_NULL;
			Mat A = PETSC_NULL;
			int num_elements, num_local_elements = local.size();
			std::vector<int> local_id(local.size());
			std::vector<double> local_phi(local.size()), rhs(local.size());
			moab::Tag global_id;
			moab::Range::iterator obj;
			moab::MeshTopoUtil topo(&imoab);
			//get local id, create phi tag, get total number of elements
			result = imoab.tag_get_handle(PARALLEL_GID_TAG_NAME,1,moab::MB_TYPE_INTEGER,global_id);
			CHKERR("cannot get global id tag\n");
			result = imoab.tag_get_data(phi_tag,local,&local_phi[0]);
			CHKERR("cannot get phi values\n");
			result = imoab.tag_get_data(global_id,local,&local_id[0]);
			CHKERR("cannot get global id values\n");
			result = MPI_Allreduce(&num_local_elements,&num_elements,1,MPI_INTEGER,MPI_SUM,MPI_COMM_WORLD);
			CHKERR("cannot do allreduce operation\n");
			start_serial_output(pcomm.comm());
			std::cout << "rank " << rank << " number of local elements " << num_local_elements << " number of global elements " << num_elements << std::endl;
			end_serial_output(pcomm.comm());
			result = PetscInitialize(&argc,&argv,PETSC_NULL,PETSC_NULL);
			CHKERR("cannot initialize petsc\n");
			result = VecCreate(PETSC_COMM_WORLD,&x);
			CHKERR("cannot create vector\n");
			result = VecSetType(x,size==1?VECSEQ:VECMPI);
			CHKERR("cannot set vector type\n");
			result = VecSetSizes(x,num_local_elements,num_elements);
			CHKERR("cannot set vector size\n");
			result = VecCreate(PETSC_COMM_WORLD,&b);
			CHKERR("cannot create vector\n");
			result = VecSetType(b,size==1?VECSEQ:VECMPI);
			CHKERR("cannot set vector type\n");
			result = VecSetSizes(b,num_local_elements,num_elements);
			CHKERR("cannot set vector size\n");
			result = MatCreate(PETSC_COMM_WORLD,&A);
			CHKERR("cannot create matrix\n");
			result = MatSetType(A,size==1?MATSEQAIJ:MATMPIAIJ);
			CHKERR("cannot set mpiaij type to matrix\n");
			result = MatSetSizes(A,num_local_elements,num_local_elements,num_elements,num_elements);
			CHKERR("cannot set size to matrix\n");
			result = MatMPIAIJSetPreallocation(A,7,PETSC_NULL,7,PETSC_NULL);
			CHKERR("cannot set preallocation for matrix\n");
			//assemble matrix
			for(its = 0, obj = local.begin(); obj != local.end(); obj++, its++)
			{
				int pos[2];
				double coords_obj[3],  obj_volume,val[2];
				moab::Range::iterator face;
				moab::Range::iterator neighbour;
				moab::Range faces;
				rhs[its] = 0.0;
				pos[0] = local_id[its];
				//result = topo.get_bridge_adjacencies(*obj,2,3,adj);
				//CHKERR("cannot get bridge adjacencies\n");
				result = topo.get_average_position(&*obj,1,coords_obj);
				CHKERR("cannot get position\n");
				result = imoab.get_adjacencies(&*obj,1,2,true,faces,moab::Interface::UNION);
				CHKERR("cannot get faces\n");
				for(face = faces.begin(); face != faces.end(); face++)
				{
					moab::Range neighbours;
					double coords_neighbour[3],face_area,h;
					result = imoab.get_adjacencies(&*face,1,3,true,neighbours,moab::Interface::UNION);
					CHKERR("cannot get face neighbours\n");
					neighbour = neighbours.begin();
					switch(neighbours.size())
					{
						case 1: //This is boundary face
							result = topo.get_average_position(&*face,1,coords_neighbour);
							CHKERR("cannot get face position\n");
							h = dist(coords_neighbour,coords_obj);
							obj_volume = h*h*h*8; //replace later with fair volume
							face_area = h*h*4; // replace later with fair area
							val[0] = -1.0 / h * face_area / obj_volume;
							if( dimdir(coords_obj,coords_neighbour) == 2 )
								rhs[its] -= 1.0 / h * face_area / obj_volume; 
							result = MatSetValues(A,1,&local_id[its],1,pos,val,ADD_VALUES);
							CHKERR("cannot set values to matrix\n");
							break;
						case 2: //This is inner face
							if( *obj == neighbours[0] )
								neighbour++;
							result = imoab.tag_get_data(global_id,&*neighbour,1,&pos[1]);
							CHKERR("cannot get id tag\n");
							result = topo.get_average_position(&*neighbour,1,coords_neighbour);
							CHKERR("cannot get position\n");
							h = dist(coords_neighbour,coords_obj);
							obj_volume = h*h*h; //replace later with fair volume
							face_area = h*h; // replace later with fair area
							val[0] = -1.0 / h * face_area / obj_volume;
							val[1] = 1.0 / h * face_area / obj_volume;
							result = MatSetValues(A,1,&local_id[its],2,pos,val,ADD_VALUES);
							CHKERR("cannot set values to matrix\n");
							break;
						default:
							std::cout << rank << " number of neighbours " << neighbours.size() << " on line " << __LINE__ << std::endl;
					}
				}
				
				
				//identity matrix
				//double one = 1.0;
				//result = MatSetValues(A,1,&local_id[its],1,&local_id[its],&one,INSERT_VALUES);
				//CHKERR("cannot set values to matrix\n");
			}
			result = VecSetValues(x,num_local_elements,&local_id[0],&local_phi[0],INSERT_VALUES);
			CHKERR("cannot set values to x vector\n");
			//for(its = 0, obj = local.begin(); obj != local.end(); obj++, its++) local_phi[its] = rand()*1.0/RAND_MAX;
			result = VecSetValues(b,num_local_elements,&local_id[0],&rhs[0],INSERT_VALUES);
			CHKERR("cannot set values to b vector\n");
			//prepare matrix and vectors
			result = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
			CHKERR("matrix assembly begin failed\n");
			result = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
			CHKERR("matrix assembly end failed\n");
			result = VecAssemblyBegin(x);
			CHKERR("vector x assembly begin failed\n");
			result = VecAssemblyEnd(x);
			CHKERR("vector x assembly end failed\n");
			result = VecAssemblyBegin(b);
			CHKERR("vector b assembly begin failed\n");
			result = VecAssemblyEnd(b);
			CHKERR("vector b assembly end failed\n");
			result = KSPCreate(PETSC_COMM_WORLD,&ksp);
			CHKERR("cannot create solver\n");
			result = KSPSetFromOptions(ksp);
			CHKERR("cannot set solver options\n");
			result = KSPSetOperators(ksp,A,A,DIFFERENT_NONZERO_PATTERN);
			CHKERR("cannot set solver operators\n");
			result = KSPSolve(ksp,b,x);
			CHKERR("solver failed\n");
			//print convergence info
			result = KSPGetIterationNumber(ksp,&its);
			CHKERR("cannot get iteration number\n");
			result = KSPGetResidualNorm(ksp,&norm);
			CHKERR("cannot get residual norm\n");
			if( rank == 0 ) std::cout << "method converged in " << its << " iterations to " << norm << " norm" << std::endl;
			//retrive info back to mesh
			result = VecGetValues(x,num_local_elements,&local_id[0],&local_phi[0]);
			CHKERR("cannot retrive data from petsc\n");
			result = imoab.tag_set_data(phi_tag,local,&local_phi[0]);
			CHKERR("cannot put data to mesh\n");
			moab::Range new_exch_tmp;
			result = pcomm.exchange_tags(phi_tag,new_exch_tmp);
			CHKERR("cannot exchange tag data\n");
			//destroy allocated vectors, matrix and solver
			result = VecDestroy(&x);
			CHKERR("cannot destroy vector x\n");
			result = VecDestroy(&b);
			CHKERR("cannot destroy vector b\n");
			result = MatDestroy(&A);
			CHKERR("cannot destroy matrix A\n");
			result = KSPDestroy(&ksp);
			CHKERR("cannot destroy solver\n");
			//finalize
			result = PetscFinalize();
			CHKERR("cannot finalize petsc\n");
		}
		
		
		{
			PIFile_write_pvtk(imoab,imoab.get_root_set(),pcomm,"out");
		}
		*/
		
	}
	else
	{
		//print help
		if( !rank ) std::cout << argv[0] << " [mesh_file]" << std::endl;
	}
	
	MPI_Finalize();
	return 0;
}