Context Navigation

source: XIOS/trunk/src/node/mesh.cpp @ 1542

Last change on this file since 1542 was 1542, checked in by oabramkina, 6 years ago

Replacing Boost's unordered_map and shared_pointer by its STL counterparts.

Two notes for Curie:

one can see the content of unordered_map with ddt only if XIOS has been compiled with gnu
XIOS will not compile any more with pgi (all available versions use old STL which are not up to the c++11 norms)

Property svn:executable set to *

File size: 82.6 KB

Rev	Line
[879]	1	/*!
	2	\file mesh.cpp
	3	\author Olga Abramkina
	4	\brief Definition of class CMesh.
	5	*/
	6
	7	#include "mesh.hpp"
[1542]	8	#include <boost/functional/hash.hpp>
	9	//#include <unordered_map>
[879]	10
	11	namespace xios {
	12
	13	/// ////////////////////// DÃ©finitions ////////////////////// ///
	14
[1158]	15	CMesh::CMesh(void) : nbNodesGlo(0), nbEdgesGlo(0)
	16	, node_start(0), node_count(0)
	17	, edge_start(0), edge_count(0)
	18	, nbFaces_(0), nbNodes_(0), nbEdges_(0)
	19	, nodesAreWritten(false), edgesAreWritten(false), facesAreWritten(false)
[879]	20	, node_lon(), node_lat()
	21	, edge_lon(), edge_lat(), edge_nodes()
	22	, face_lon(), face_lat()
	23	, face_nodes()
[1158]	24	, pNodeGlobalIndex(NULL), pEdgeGlobalIndex(NULL)
[879]	25	{
	26	}
	27
	28
	29	CMesh::~CMesh(void)
	30	{
[924]	31	if (pNodeGlobalIndex != NULL) delete pNodeGlobalIndex;
	32	if (pEdgeGlobalIndex != NULL) delete pEdgeGlobalIndex;
[879]	33	}
	34
[881]	35	std::map <StdString, CMesh> CMesh::meshList = std::map <StdString, CMesh>();
[924]	36	std::map <StdString, vector<int> > CMesh::domainList = std::map <StdString, vector<int> >();
[879]	37
	38	///---------------------------------------------------------------
	39	/*!
[881]	40	* \fn bool CMesh::getMesh (StdString meshName)
	41	* Returns a pointer to a mesh. If a mesh has not been created, creates it and adds its name to the list of meshes meshList.
[879]	42	* \param [in] meshName The name of a mesh ("name" attribute of a domain).
[924]	43	* \param [in] nvertex Number of verteces (1 for nodes, 2 for edges, 3 and up for faces).
[879]	44	*/
[924]	45	CMesh* CMesh::getMesh (StdString meshName, int nvertex)
[879]	46	{
[924]	47	CMesh::domainList[meshName].push_back(nvertex);
	48
[883]	49	if ( CMesh::meshList.begin() != CMesh::meshList.end() )
[879]	50	{
[883]	51	for (std::map<StdString, CMesh>::iterator it=CMesh::meshList.begin(); it!=CMesh::meshList.end(); ++it)
	52	{
	53	if (it->first == meshName)
	54	return &meshList[meshName];
	55	else
	56	{
	57	CMesh newMesh;
	58	CMesh::meshList.insert( make_pair(meshName, newMesh) );
	59	return &meshList[meshName];
	60	}
	61	}
	62	}
	63	else
	64	{
[882]	65	CMesh newMesh;
[881]	66	CMesh::meshList.insert( make_pair(meshName, newMesh) );
[883]	67	return &meshList[meshName];
[879]	68	}
	69	}
	70
	71	///----------------------------------------------------------------
[924]	72	size_t hashPair(size_t first, size_t second)
[879]	73	{
[924]	74	HashXIOS<size_t> sizetHash;
	75	size_t seed = sizetHash(first) + 0x9e3779b9 ;
	76	seed ^= sizetHash(second) + 0x9e3779b9 + (seed << 6) + (seed >> 2);
[879]	77	return seed ;
	78	}
	79
	80	///----------------------------------------------------------------
[924]	81	size_t hashPairOrdered(size_t first, size_t second)
	82	{
	83	size_t seed;
	84	HashXIOS<size_t> sizetHash;
	85	if (first < second)
	86	{
	87	seed = sizetHash(first) + 0x9e3779b9 ;
	88	seed ^= sizetHash(second) + 0x9e3779b9 + (seed << 6) + (seed >> 2);
	89	}
	90	else
	91	{
	92	seed = sizetHash(second) + 0x9e3779b9 ;
	93	seed ^= sizetHash(first) + 0x9e3779b9 + (seed << 6) + (seed >> 2);
	94	}
	95	return seed ;
	96	}
	97
	98	///----------------------------------------------------------------
[879]	99	/*!
[924]	100	* \fn size_t generateNodeIndex(vector<size_t>& valList, int rank)
	101	* Generates a node index.
	102	* If the same node is generated by two processes, each process will have its own node index.
	103	* \param [in] valList Vector storing four node hashes.
	104	* \param [in] rank MPI process rank.
	105	*/
	106	size_t generateNodeIndex(vector<size_t>& valList, int rank)
	107	{
	108	// Sort is needed to avoid problems for nodes with lon = 0 generated by faces in east and west semisphere
	109	vector<size_t> vec = valList;
	110	sort (vec.begin(), vec.end());
	111	size_t seed = rank ;
	112	int it = 0;
	113	for(; it != vec.size(); ++it)
	114	{
	115	seed = hashPair(seed, vec[it]);
	116	}
	117	return seed ;
	118	}
	119
[1158]	120	///----------------------------------------------------------------
	121	/*!
	122	* \fn size_t generateNodeIndex(vector<size_t>& valList)
	123	* Generates a node index unique for all processes.
	124	* \param [in] valList Vector storing four node hashes.
	125	*/
	126	size_t generateNodeIndex(vector<size_t>& valList)
	127	{
	128	// Sort is needed to avoid problems for nodes with lon = 0 generated by faces in east and west semisphere
	129	vector<size_t> vec = valList;
	130	sort (vec.begin(), vec.end());
	131	size_t seed = vec[0] ;
	132	int it = 1;
	133	for(; it != vec.size(); ++it)
	134	{
	135	seed = hashPair(seed, vec[it]);
	136	}
	137	return seed ;
	138	}
	139
[924]	140	///----------------------------------------------------------------
	141	/*!
	142	* \fn CArray<size_t,1>& CMesh::createHashes (const double longitude, const double latitude)
[879]	143	* Creates two hash values for each dimension, longitude and latitude.
[924]	144	* \param [in] longitude Node longitude in degrees.
	145	* \param [in] latitude Node latitude in degrees ranged from 0 to 360.
[879]	146	*/
[881]	147
[924]	148	vector<size_t> CMesh::createHashes (const double longitude, const double latitude)
[879]	149	{
	150	double minBoundLon = 0. ;
	151	double maxBoundLon = 360. ;
[924]	152	double minBoundLat = -90. ;
	153	double maxBoundLat = 90. ;
[879]	154	double prec=1e-11 ;
	155	double precLon=prec ;
	156	double precLat=prec ;
[924]	157	double lon = longitude;
	158	double lat = latitude;
[879]	159
[924]	160	if (lon > (360.- prec)) lon = 0.;
	161
[879]	162	size_t maxsize_t=numeric_limits<size_t>::max() ;
	163	if ( (maxBoundLon-minBoundLon)/maxsize_t > precLon) precLon=(maxBoundLon-minBoundLon)/maxsize_t ;
	164	if ( (maxBoundLat-minBoundLat)/maxsize_t > precLat) precLat=(maxBoundLat-minBoundLat)/maxsize_t ;
	165
	166	size_t iMinLon=0 ;
	167	size_t iMaxLon=(maxBoundLon-minBoundLon)/precLon ;
	168	size_t iMinLat=0 ;
	169	size_t iMaxLat=(maxBoundLat-minBoundLat)/precLat ;
	170
[900]	171	vector<size_t> hash(4);
[879]	172	size_t lon0,lon1,lat0,lat1 ;
	173
	174	lon0=(lon-minBoundLon)/precLon ;
	175	if ( ((lon0+1)precLon + lon0precLon)/2 > lon-minBoundLon)
	176	{
	177	if (lon0==iMinLon) lon1=iMaxLon ;
	178	else lon1=lon0-1 ;
	179	}
	180	else
	181	{
	182	if (lon0==iMaxLon) lon1=iMinLon ;
	183	else lon1=lon0+1 ;
	184	}
	185
	186	lat0=(lat-minBoundLat)/precLat ;
	187	if ( ((lat0+1)precLat + lat0precLat)/2 > lat-minBoundLat)
	188	{
	189	if (lat0==iMinLat) lat1=lat0 ;
	190	else lat1=lat0-1 ;
	191	}
	192	else
	193	{
	194	if (lat0==iMaxLat) lat1=lat0 ;
	195	else lat1=lat0+1 ;
	196	}
	197
[900]	198	hash[0] = hashPair(lon0,lat0) ;
	199	hash[1] = hashPair(lon0,lat1) ;
	200	hash[2] = hashPair(lon1,lat0) ;
	201	hash[3] = hashPair(lon1,lat1) ;
[879]	202
[900]	203	return hash;
[879]	204
[900]	205	} // createHashes
	206
[879]	207	///----------------------------------------------------------------
	208	std::pair<int,int> make_ordered_pair(int a, int b)
	209	{
	210	if ( a < b )
	211	return std::pair<int,int>(a,b);
	212	else
	213	return std::pair<int,int>(b,a);
	214	}
	215
	216	///----------------------------------------------------------------
	217	/*!
	218	* \fn void CMesh::createMesh(const CArray<double, 1>& lonvalue, const CArray<double, 1>& latvalue,
	219	const CArray<double, 2>& bounds_lon, const CArray<double, 2>& bounds_lat)
	220	* Creates or updates a mesh for the three types of mesh elements: nodes, edges, and faces.
	221	* \param [in] lonvalue Array of longitudes.
	222	* \param [in] latvalue Array of latitudes.
	223	* \param [in] bounds_lon Array of boundary longitudes. Its size depends on the element type.
[881]	224	* \param [in] bounds_lat Array of boundary latitudes. Its size depends on the element type.
[879]	225	*/
[1542]	226	// void CMesh::createMesh(const CArray<double, 1>& lonvalue, const CArray<double, 1>& latvalue,
	227	// const CArray<double, 2>& bounds_lon, const CArray<double, 2>& bounds_lat)
	228	// {
	229	// int nvertex = (bounds_lon.numElements() == 0) ? 1 : bounds_lon.rows();
	230	//
	231	// if (nvertex == 1)
	232	// {
	233	// nbNodes_ = lonvalue.numElements();
	234	// node_lon.resizeAndPreserve(nbNodes_);
	235	// node_lat.resizeAndPreserve(nbNodes_);
	236	// for (int nn = 0; nn < nbNodes_; ++nn)
	237	// {
	238	// if (map_nodes.find(make_pair (lonvalue(nn), latvalue(nn))) == map_nodes.end())
	239	// {
	240	// map_nodes[make_pair (lonvalue(nn), latvalue(nn))] = nn ;
	241	// node_lon(nn) = lonvalue(nn);
	242	// node_lat(nn) = latvalue(nn);
	243	// }
	244	// }
	245	// }
	246	// else if (nvertex == 2)
	247	// {
	248	// nbEdges_ = bounds_lon.shape()[1];
	249	//
	250	// // Create nodes and edge_node connectivity
	251	// node_lon.resizeAndPreserve(nbEdges_*nvertex); // Max possible number of nodes
	252	// node_lat.resizeAndPreserve(nbEdges_*nvertex);
	253	// edge_nodes.resizeAndPreserve(nvertex, nbEdges_);
	254	//
	255	// for (int ne = 0; ne < nbEdges_; ++ne)
	256	// {
	257	// for (int nv = 0; nv < nvertex; ++nv)
	258	// {
	259	// if (map_nodes.find(make_pair (bounds_lon(nv, ne), bounds_lat(nv ,ne))) == map_nodes.end())
	260	// {
	261	// map_nodes[make_pair (bounds_lon(nv, ne), bounds_lat(nv, ne))] = nbNodes_ ;
	262	// edge_nodes(nv,ne) = nbNodes_ ;
	263	// node_lon(nbNodes_) = bounds_lon(nv, ne);
	264	// node_lat(nbNodes_) = bounds_lat(nv, ne);
	265	// ++nbNodes_ ;
	266	// }
	267	// else
	268	// edge_nodes(nv,ne) = map_nodes[make_pair (bounds_lon(nv, ne), bounds_lat(nv ,ne))];
	269	// }
	270	// }
	271	// node_lon.resizeAndPreserve(nbNodes_);
	272	// node_lat.resizeAndPreserve(nbNodes_);
	273	//
	274	// // Create edges
	275	// edge_lon.resizeAndPreserve(nbEdges_);
	276	// edge_lat.resizeAndPreserve(nbEdges_);
	277	//
	278	// for (int ne = 0; ne < nbEdges_; ++ne)
	279	// {
	280	// if (map_edges.find(make_ordered_pair (edge_nodes(0,ne), edge_nodes(1,ne))) == map_edges.end())
	281	// {
	282	// map_edges[make_ordered_pair ( edge_nodes(0,ne), edge_nodes(1,ne) )] = ne ;
	283	// edge_lon(ne) = lonvalue(ne);
	284	// edge_lat(ne) = latvalue(ne);
	285	// }
	286	//
	287	// }
	288	// edgesAreWritten = true;
	289	// }
	290	// else
	291	// {
	292	// nbFaces_ = bounds_lon.shape()[1];
	293	//
	294	// // Create nodes and face_node connectivity
	295	// node_lon.resizeAndPreserve(nbFaces_*nvertex); // Max possible number of nodes
	296	// node_lat.resizeAndPreserve(nbFaces_*nvertex);
	297	// face_nodes.resize(nvertex, nbFaces_);
	298	//
	299	// for (int nf = 0; nf < nbFaces_; ++nf)
	300	// {
	301	// for (int nv = 0; nv < nvertex; ++nv)
	302	// {
	303	// if (map_nodes.find(make_pair (bounds_lon(nv, nf), bounds_lat(nv ,nf))) == map_nodes.end())
	304	// {
	305	// map_nodes[make_pair (bounds_lon(nv, nf), bounds_lat(nv, nf))] = nbNodes_ ;
	306	// face_nodes(nv,nf) = nbNodes_ ;
	307	// node_lon(nbNodes_) = bounds_lon(nv, nf);
	308	// node_lat(nbNodes_) = bounds_lat(nv ,nf);
	309	// ++nbNodes_ ;
	310	// }
	311	// else
	312	// {
	313	// face_nodes(nv,nf) = map_nodes[make_pair (bounds_lon(nv, nf), bounds_lat(nv ,nf))];
	314	// }
	315	// }
	316	// }
	317	// node_lon.resizeAndPreserve(nbNodes_);
	318	// node_lat.resizeAndPreserve(nbNodes_);
	319	//
	320	// // Create edges and edge_nodes connectivity
	321	// edge_lon.resizeAndPreserve(nbFaces_*nvertex); // Max possible number of edges
	322	// edge_lat.resizeAndPreserve(nbFaces_*nvertex);
	323	// edge_nodes.resizeAndPreserve(2, nbFaces_*nvertex);
	324	// edge_faces.resize(2, nbFaces_*nvertex);
	325	// face_edges.resize(nvertex, nbFaces_);
	326	// face_faces.resize(nvertex, nbFaces_);
	327	//
	328	// vector<int> countEdges(nbFaces_*nvertex); // needed in case if edges have been already generated
	329	// vector<int> countFaces(nbFaces_);
	330	// countEdges.assign(nbFaces_*nvertex, 0);
	331	// countFaces.assign(nbFaces_, 0);
	332	// int edge;
	333	// for (int nf = 0; nf < nbFaces_; ++nf)
	334	// {
	335	// for (int nv1 = 0; nv1 < nvertex; ++nv1)
	336	// {
	337	// int nv = 0;
	338	// int nv2 = (nv1 < nvertex -1 ) ? (nv1 + 1) : (nv1 + 1 - nvertex); // cyclic rotation
	339	// if (map_edges.find(make_ordered_pair (face_nodes(nv1,nf), face_nodes(nv2,nf))) == map_edges.end())
	340	// {
	341	// map_edges[make_ordered_pair (face_nodes(nv1,nf), face_nodes(nv2,nf))] = nbEdges_ ;
	342	// face_edges(nv1,nf) = map_edges[make_ordered_pair (face_nodes(nv1,nf), face_nodes(nv2,nf))];
	343	// edge_faces(0,nbEdges_) = nf;
	344	// edge_faces(1,nbEdges_) = -999;
	345	// face_faces(nv1,nf) = 999999;
	346	// edge_nodes(Range::all(),nbEdges_) = face_nodes(nv1,nf), face_nodes(nv2,nf);
	347	// edge_lon(nbEdges_) = ( abs( node_lon(face_nodes(nv1,nf)) - node_lon(face_nodes(nv2,nf))) < 180.) ?
	348	// (( node_lon(face_nodes(nv1,nf)) + node_lon(face_nodes(nv2,nf))) * 0.5) :
	349	// (( node_lon(face_nodes(nv1,nf)) + node_lon(face_nodes(nv2,nf))) * 0.5 -180.);
	350	// edge_lat(nbEdges_) = ( node_lat(face_nodes(nv1,nf)) + node_lat(face_nodes(nv2,nf)) ) * 0.5;
	351	// ++nbEdges_;
	352	// }
	353	// else
	354	// {
	355	// edge = map_edges[make_ordered_pair (face_nodes(nv1,nf), face_nodes(nv2,nf))];
	356	// face_edges(nv1,nf) = edge;
	357	// if (edgesAreWritten)
	358	// {
	359	// edge_faces(countEdges[edge], edge) = nf;
	360	// if (countEdges[edge]==0)
	361	// {
	362	// face_faces(nv1,nf) = 999999;
	363	// }
	364	// else
	365	// {
	366	// int face1 = nf; // = edge_faces(1,edge)
	367	// int face2 = edge_faces(0,edge);
	368	// face_faces(countFaces[face1], face1) = face2;
	369	// face_faces(countFaces[face2], face2) = face1;
	370	// ++(countFaces[face1]);
	371	// ++(countFaces[face2]);
	372	// }
	373	// }
	374	// else
	375	// {
	376	// edge_faces(1,edge) = nf;
	377	// int face1 = nf; // = edge_faces(1,edge)
	378	// int face2 = edge_faces(0,edge);
	379	// face_faces(countFaces[face1], face1) = face2;
	380	// face_faces(countFaces[face2], face2) = face1;
	381	// ++(countFaces[face1]);
	382	// ++(countFaces[face2]);
	383	// }
	384	// ++(countEdges[edge]);
	385	// }
	386	// }
	387	// }
	388	// edge_nodes.resizeAndPreserve(2, nbEdges_);
	389	// edge_faces.resizeAndPreserve(2, nbEdges_);
	390	// edge_lon.resizeAndPreserve(nbEdges_);
	391	// edge_lat.resizeAndPreserve(nbEdges_);
	392	//
	393	// // Create faces
	394	// face_lon.resize(nbFaces_);
	395	// face_lat.resize(nbFaces_);
	396	// face_lon = lonvalue;
	397	// face_lat = latvalue;
	398	// facesAreWritten = true;
	399	//
	400	// } // nvertex > 2
	401	//
	402	// } // createMesh()
[879]	403
	404	///----------------------------------------------------------------
	405	/*!
	406	* \fn void CMesh::createMeshEpsilon(const CArray<double, 1>& lonvalue, const CArray<double, 1>& latvalue,
	407	const CArray<double, 2>& bounds_lon, const CArray<double, 2>& bounds_lat)
	408	* Creates or updates a mesh for the three types of mesh elements: nodes, edges, and faces.
	409	* Precision check is implemented with two hash values for each dimension, longitude and latitude.
[924]	410	* \param [in] comm
[879]	411	* \param [in] lonvalue Array of longitudes.
	412	* \param [in] latvalue Array of latitudes.
	413	* \param [in] bounds_lon Array of boundary longitudes. Its size depends on the element type.
[881]	414	* \param [in] bounds_lat Array of boundary latitudes. Its size depends on the element type.
[879]	415	*/
[924]	416	void CMesh::createMeshEpsilon(const MPI_Comm& comm,
	417	const CArray<double, 1>& lonvalue, const CArray<double, 1>& latvalue,
	418	const CArray<double, 2>& bounds_lon, const CArray<double, 2>& bounds_lat)
[879]	419	{
[900]	420
[924]	421	int nvertex = (bounds_lon.numElements() == 0) ? 1 : bounds_lon.rows();
	422	int mpiRank, mpiSize;
	423	MPI_Comm_rank(comm, &mpiRank);
	424	MPI_Comm_size(comm, &mpiSize);
[929]	425	double prec = 1e-11; // used in calculations of edge_lon/lat
[879]	426
	427	if (nvertex == 1)
	428	{
[929]	429	nbNodes_ = lonvalue.numElements();
	430	node_lon.resize(nbNodes_);
	431	node_lat.resize(nbNodes_);
[924]	432	node_lon = lonvalue;
	433	node_lat = latvalue;
[900]	434
[924]	435	// Global node indexes
	436	vector<size_t> hashValues(4);
	437	CClientClientDHTSizet::Index2VectorInfoTypeMap nodeHash2IdxGlo;
[929]	438	for (size_t nn = 0; nn < nbNodes_; ++nn)
[879]	439	{
[900]	440	hashValues = CMesh::createHashes(lonvalue(nn), latvalue(nn));
[924]	441	for (size_t nh = 0; nh < 4; ++nh)
	442	{
[929]	443	nodeHash2IdxGlo[hashValues[nh]].push_back(mpiRank*nbNodes_ + nn);
[924]	444	}
[879]	445	}
[924]	446	pNodeGlobalIndex = new CClientClientDHTSizet (nodeHash2IdxGlo, comm);
	447	nodesAreWritten = true;
	448	}
[900]	449
[924]	450	else if (nvertex == 2)
	451	{
[929]	452	nbEdges_ = bounds_lon.shape()[1];
	453	edge_lon.resize(nbEdges_);
	454	edge_lat.resize(nbEdges_);
[924]	455	edge_lon = lonvalue;
	456	edge_lat = latvalue;
[929]	457	edge_nodes.resize(nvertex, nbEdges_);
	458
	459	// For determining the global edge index
[1507]	460	unsigned long nbEdgesOnProc = nbEdges_;
	461	unsigned long nbEdgesAccum;
[929]	462	MPI_Scan(&nbEdgesOnProc, &nbEdgesAccum, 1, MPI_UNSIGNED_LONG, MPI_SUM, comm);
	463	nbEdgesAccum -= nbEdges_;
	464
[924]	465	CClientClientDHTSizet::Index2VectorInfoTypeMap edgeHash2IdxGlo;
	466	CClientClientDHTSizet::Index2VectorInfoTypeMap edgeHash2Idx;
[900]	467
[924]	468	// Case (1): node indexes have been generated by domain "nodes"
	469	if (nodesAreWritten)
	470	{
	471	vector<size_t> hashValues(4);
[929]	472	CArray<size_t,1> nodeHashList(nbEdges_nvertex4);
	473	for (int ne = 0; ne < nbEdges_; ++ne) // size_t doesn't work with CArray<double, 2>
[924]	474	{
	475	for (int nv = 0; nv < nvertex; ++nv)
	476	{
	477	hashValues = CMesh::createHashes(bounds_lon(nv, ne), bounds_lat(nv, ne));
	478	for (int nh = 0; nh < 4; ++nh)
	479	{
	480	nodeHashList((nenvertex + nv)4 + nh) = hashValues[nh];
	481	}
	482	}
	483	}
[900]	484
[924]	485	// Recuperating the node global indexing and writing edge_nodes
	486	// Creating map edgeHash2IdxGlo <hash, idxGlo>
	487	pNodeGlobalIndex->computeIndexInfoMapping(nodeHashList);
	488	CClientClientDHTSizet::Index2VectorInfoTypeMap& nodeHash2IdxGlo = pNodeGlobalIndex->getInfoIndexMap();
	489	CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator it;
[929]	490	size_t nodeIdxGlo1, nodeIdxGlo2;
	491	for (int ne = 0; ne < nbEdges_; ++ne)
[924]	492	{
	493	for (int nv = 0; nv < nvertex; ++nv)
	494	{
	495	int nh = 0;
	496	it = nodeHash2IdxGlo.find(nodeHashList((nenvertex + nv)4 + nh));
	497	// The loop below is needed in case if a hash generated by domain "edges" differs
	498	// from that generated by domain "nodes" because of possible precision issues
	499	while (it == nodeHash2IdxGlo.end())
	500	{
	501	++nh;
	502	it = nodeHash2IdxGlo.find(nodeHashList((nenvertex + nv)4 + nh));
	503	}
	504	edge_nodes(nv,ne) = it->second[0];
[929]	505	if (nv ==0)
	506	nodeIdxGlo1 = it->second[0];
	507	else
	508	nodeIdxGlo2 = it->second[0];
[924]	509	}
[929]	510	size_t edgeIdxGlo = nbEdgesAccum + ne;
	511	edgeHash2IdxGlo[ hashPairOrdered(nodeIdxGlo1, nodeIdxGlo2) ].push_back(edgeIdxGlo);
[924]	512	}
	513	} // nodesAreWritten
[900]	514
[929]	515
[924]	516	// Case (2): node indexes have not been generated previously
	517	else
	518	{
	519	// (2.1) Creating a list of hashes for each node and a map nodeHash2Idx <hash, <idx,rank> >
	520	vector<size_t> hashValues(4);
	521	CClientClientDHTSizet::Index2VectorInfoTypeMap nodeHash2Idx;
[929]	522	CArray<size_t,1> nodeHashList(nbEdges_nvertex4);
[1158]	523	int nbHash = 0;
[929]	524	for (int ne = 0; ne < nbEdges_; ++ne)
[924]	525	{
	526	for (int nv = 0; nv < nvertex; ++nv)
	527	{
	528	hashValues = CMesh::createHashes(bounds_lon(nv, ne), bounds_lat(nv, ne));
	529	for (int nh = 0; nh < 4; ++nh)
	530	{
	531	if (nodeHash2Idx[hashValues[nh]].size() == 0)
	532	{
[1158]	533	nodeHash2Idx[hashValues[nh]].push_back(generateNodeIndex(hashValues));
[924]	534	nodeHash2Idx[hashValues[nh]].push_back(mpiRank);
[1158]	535	nodeHashList(nbHash) = hashValues[nh];
	536	++nbHash;
[924]	537	}
	538	}
	539	}
	540	}
[1158]	541	nodeHashList.resizeAndPreserve(nbHash);
[900]	542
[924]	543	// (2.2) Generating global node indexes
[1158]	544	// The ownership criterion: priority of the process of smaller index
[924]	545	// Maps generated in this step are:
[1158]	546	// Maps generated in this step are:
	547	// nodeHash2Info = <hash, [[idx, rankMin], [idx, rank1], [idx, rank3]..]>
	548	// nodeIdx2Idx = <idx, <rankOwner, idx>>
[900]	549
[1158]	550	CClientClientDHTSizet dhtNodeHash(nodeHash2Idx, comm);
	551	dhtNodeHash.computeIndexInfoMapping(nodeHashList);
	552	CClientClientDHTSizet::Index2VectorInfoTypeMap& nodeHash2Info = dhtNodeHash.getInfoIndexMap();
[900]	553
	554
[1158]	555	CClientClientDHTSizet::Index2VectorInfoTypeMap nodeIdx2Idx;
	556	CArray<size_t,1> nodeIdxList(nbEdges_nvertex4);
	557	size_t nIdx = 0;
[900]	558
[1158]	559	for (CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator it = nodeHash2Info.begin(); it != nodeHash2Info.end(); ++it)
	560	{
	561	size_t rankMin = (it->second)[1];
	562	size_t idx = (it->second)[0];
	563	for (int i = 2; i < (it->second).size();)
	564	{
	565	if ( (it->second)[i+1] < rankMin)
	566	{
	567	idx = (it->second)[i];
	568	rankMin = (it->second)[i+1];
	569	(it->second)[i+1] = (it->second)[i-1];
	570	}
	571	i += 2;
	572	}
	573	if (nodeIdx2Idx.count(idx) == 0)
	574	{
	575	if (mpiRank == rankMin)
	576	{
	577	nodeIdx2Idx[idx].push_back(rankMin);
	578	nodeIdx2Idx[idx].push_back(idx);
	579	}
	580	nodeIdxList(nIdx) = idx;
	581	++nIdx;
	582	}
	583	}
	584	nodeIdxList.resizeAndPreserve(nIdx);
	585
	586	// CDHTAutoIndexing will not give consistent node numbering for varying number of procs. =>
	587	// Solution: global node indexing by hand.
	588	// Maps modified in this step:
	589	// nodeIdx2Idx = <idx, idxGlo>
[1507]	590	unsigned long nodeCount = nodeIdx2Idx.size();
	591	unsigned long nodeStart, nbNodes;
[1158]	592	MPI_Scan(&nodeCount, &nodeStart, 1, MPI_UNSIGNED_LONG, MPI_SUM, comm);
	593	int nNodes = nodeStart;
	594	MPI_Bcast(&nNodes, 1, MPI_UNSIGNED_LONG, mpiSize-1, comm);
	595	nbNodesGlo = nNodes;
	596
	597	nodeStart -= nodeCount;
	598	node_start = nodeStart;
	599	node_count = nodeCount;
	600	CClientClientDHTSizet::Index2VectorInfoTypeMap dummyMap; // just a dummy map used to ensure that each node is numbered only once
	601	size_t count = 0;
	602
	603	for (int ne = 0; ne < nbEdges_; ++ne)
	604	{
	605	for (int nv = 0; nv < nvertex; ++nv)
	606	{
	607	vector<size_t> hashValues = CMesh::createHashes(bounds_lon(nv, ne), bounds_lat(nv, ne));
	608	size_t nodeIdx = generateNodeIndex(hashValues);
	609	CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator it = nodeIdx2Idx.find(nodeIdx);
	610	if (it != nodeIdx2Idx.end())
	611	{
	612	if (dummyMap.count(nodeIdx) == 0)
	613	{
	614	dummyMap[nodeIdx].push_back(nodeIdx);
	615	(it->second)[1] = node_start + count;
	616	++count;
	617	}
	618	}
	619	}
	620	}
	621
	622	CClientClientDHTSizet dhtNodeIdx(nodeIdx2Idx, comm);
	623	dhtNodeIdx.computeIndexInfoMapping(nodeIdxList);
	624	CClientClientDHTSizet::Index2VectorInfoTypeMap& nodeIdx2IdxGlo = dhtNodeIdx.getInfoIndexMap();
	625
[924]	626	// (2.3) Saving variables: node_lon, node_lat, edge_nodes
	627	// Creating map nodeHash2IdxGlo <hash, idxGlo>
	628	// Creating map edgeHash2IdxGlo <hash, idxGlo>
[1158]	629	// nbNodesGlo = dhtNodeIdxGlo.getNbIndexesGlobal();
	630	// node_count = dhtNodeIdxGlo.getIndexCount();
	631	// node_start = dhtNodeIdxGlo.getIndexStart();
[924]	632	CClientClientDHTSizet::Index2VectorInfoTypeMap nodeHash2IdxGlo;
	633	node_lon.resize(node_count);
	634	node_lat.resize(node_count);
	635	vector <size_t> edgeNodes;
	636	size_t idxGlo = 0;
	637
[929]	638	for (int ne = 0; ne < nbEdges_; ++ne)
[924]	639	{
	640	for (int nv = 0; nv < nvertex; ++nv)
	641	{
	642	hashValues = CMesh::createHashes(bounds_lon(nv, ne), bounds_lat(nv, ne));
[1158]	643	size_t myIdx = generateNodeIndex(hashValues);
	644	CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itIdx = nodeIdx2IdxGlo.find(myIdx);
	645	idxGlo = (itIdx->second)[1];
[924]	646
[1158]	647	if (mpiRank == (itIdx->second)[0])
[924]	648	{
	649	// node_lon(idxGlo - node_start) = (bounds_lon(nv, ne) == 360.) ? (0.) : (bounds_lon(nv, ne));
	650	node_lon(idxGlo - node_start) = bounds_lon(nv, ne);
	651	node_lat(idxGlo - node_start) = bounds_lat(nv, ne);
	652	}
	653	edge_nodes(nv,ne) = idxGlo;
	654	for (int nh = 0; nh < 4; ++nh)
	655	nodeHash2IdxGlo[hashValues[nh]].push_back(idxGlo);
	656	edgeNodes.push_back(idxGlo);
	657	}
[929]	658	if (edgeNodes[0] != edgeNodes[1])
	659	{
	660	size_t edgeIdxGlo = nbEdgesAccum + ne;
	661	edgeHash2IdxGlo[ hashPairOrdered(edgeNodes[0], edgeNodes[1]) ].push_back(edgeIdxGlo);
	662	}
[924]	663	edgeNodes.clear();
	664	}
	665	pNodeGlobalIndex = new CClientClientDHTSizet (nodeHash2IdxGlo, comm);
	666	} // !nodesAreWritten
	667
	668	pEdgeGlobalIndex = new CClientClientDHTSizet (edgeHash2IdxGlo, comm);
	669	edgesAreWritten = true;
	670	} //nvertex = 2
	671
	672	else
[879]	673	{
[929]	674	nbFaces_ = bounds_lon.shape()[1];
	675	face_lon.resize(nbFaces_);
	676	face_lat.resize(nbFaces_);
[924]	677	face_lon = lonvalue;
	678	face_lat = latvalue;
[929]	679	face_nodes.resize(nvertex, nbFaces_);
	680	face_edges.resize(nvertex, nbFaces_);
[879]	681
[929]	682	// For determining the global face index
[1507]	683	unsigned long nbFacesOnProc = nbFaces_;
	684	unsigned long nbFacesAccum;
[929]	685	MPI_Scan(&nbFacesOnProc, &nbFacesAccum, 1, MPI_UNSIGNED_LONG, MPI_SUM, comm);
	686	nbFacesAccum -= nbFaces_;
	687
[924]	688	// Case (1): edges have been previously generated
	689	if (edgesAreWritten)
[879]	690	{
[924]	691	// (1.1) Recuperating node global indexing and saving face_nodes
	692	vector<size_t> hashValues(4);
[929]	693	CArray<size_t,1> nodeHashList(nbFaces_nvertex4);
	694	for (int nf = 0; nf < nbFaces_; ++nf)
[879]	695	{
[924]	696	for (int nv = 0; nv < nvertex; ++nv)
	697	{
	698	hashValues = CMesh::createHashes(bounds_lon(nv, nf), bounds_lat(nv, nf));
	699	for (int nh = 0; nh < 4; ++nh)
	700	nodeHashList((nfnvertex + nv)4 + nh) = hashValues[nh];
	701	}
	702	}
	703	pNodeGlobalIndex->computeIndexInfoMapping(nodeHashList);
	704	CClientClientDHTSizet::Index2VectorInfoTypeMap& nodeHash2IdxGlo = pNodeGlobalIndex->getInfoIndexMap();
	705	CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator it1, it2;
[929]	706	CArray<size_t,1> edgeHashList(nbFaces_*nvertex);
	707	size_t nEdge = 0;
	708	for (int nf = 0; nf < nbFaces_; ++nf)
[924]	709	{
	710	for (int nv1 = 0; nv1 < nvertex; ++nv1)
	711	{
	712	int nh1 = 0;
	713	int nv2 = (nv1 < nvertex -1 ) ? (nv1 + 1) : (nv1 + 1 - nvertex); // cyclic rotation
	714	it1 = nodeHash2IdxGlo.find(nodeHashList((nfnvertex + nv1)4 + nh1));
	715	while (it1 == nodeHash2IdxGlo.end())
	716	{
	717	++nh1;
	718	it1 = nodeHash2IdxGlo.find(nodeHashList((nfnvertex + nv1)4 + nh1));
	719	}
	720	int nh2 = 0;
	721	it2 = nodeHash2IdxGlo.find(nodeHashList((nfnvertex + nv2)4 + nh2));
	722	while (it2 == nodeHash2IdxGlo.end())
	723	{
	724	++nh2;
	725	it2 = nodeHash2IdxGlo.find(nodeHashList((nfnvertex + nv1)4 + nh2));
	726	}
	727	face_nodes(nv1,nf) = it1->second[0];
[929]	728	if (it1->second[0] != it2->second[0])
	729	{
	730	edgeHashList(nEdge) = hashPairOrdered(it1->second[0], it2->second[0]);
	731	++nEdge;
	732	}
[924]	733	}
	734	}
[929]	735	edgeHashList.resizeAndPreserve(nEdge);
[900]	736
[924]	737	// (1.2) Recuperating edge global indexing and saving face_edges
	738	pEdgeGlobalIndex->computeIndexInfoMapping(edgeHashList);
	739	CClientClientDHTSizet::Index2VectorInfoTypeMap& edgeHash2IdxGlo = pEdgeGlobalIndex->getInfoIndexMap();
	740	CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itEdgeHash;
	741	CClientClientDHTSizet::Index2VectorInfoTypeMap edgeHash2Rank;
	742	CClientClientDHTSizet::Index2VectorInfoTypeMap edgeIdxGlo2Face;
[1158]	743	CArray<size_t,1> edgeIdxList(nbFaces_*nvertex);
[924]	744	size_t iIdx = 0;
	745
[929]	746	for (int nf = 0; nf < nbFaces_; ++nf)
[924]	747	{
	748	for (int nv1 = 0; nv1 < nvertex; ++nv1)
	749	{
	750	int nh1 = 0;
	751	int nv2 = (nv1 < nvertex -1 ) ? (nv1 + 1) : (nv1 + 1 - nvertex); // cyclic rotation
	752	it1 = nodeHash2IdxGlo.find(nodeHashList((nfnvertex + nv1)4 + nh1));
	753	while (it1 == nodeHash2IdxGlo.end())
	754	{
	755	++nh1;
	756	it1 = nodeHash2IdxGlo.find(nodeHashList((nfnvertex + nv1)4 + nh1));
	757	}
	758	int nh2 = 0;
	759	it2 = nodeHash2IdxGlo.find(nodeHashList((nfnvertex + nv2)4 + nh2));
	760	while (it2 == nodeHash2IdxGlo.end())
	761	{
	762	++nh2;
	763	it2 = nodeHash2IdxGlo.find(nodeHashList((nfnvertex + nv1)4 + nh2));
	764	}
[929]	765	if (it1->second[0] != it2->second[0])
[924]	766	{
[929]	767	size_t faceIdxGlo = nbFacesAccum + nf;
	768	size_t edgeHash = hashPairOrdered(it1->second[0], it2->second[0]);
	769	itEdgeHash = edgeHash2IdxGlo.find(edgeHash);
	770	size_t edgeIdxGlo = itEdgeHash->second[0];
	771	face_edges(nv1,nf) = edgeIdxGlo;
	772	if (edgeIdxGlo2Face.count(edgeIdxGlo) == 0)
	773	{
[1158]	774	edgeIdxList(iIdx) = edgeIdxGlo;
[929]	775	++iIdx;
	776	}
	777	edgeIdxGlo2Face[edgeIdxGlo].push_back(faceIdxGlo);
	778	edgeHash2Rank[edgeHash].push_back(mpiRank);
[1158]	779	edgeHash2Rank[edgeHash].push_back(itEdgeHash->second[0]);
[924]	780	}
[929]	781	else
	782	{
	783	face_edges(nv1,nf) = 999999;
	784	}
[924]	785	}
[879]	786	}
[1158]	787	edgeIdxList.resizeAndPreserve(iIdx);
[879]	788
[924]	789	// (1.3) Saving remaining variables edge_faces and face_faces
[900]	790
[924]	791	// Establishing edge ownership
	792	// The ownership criterion: priority of the process with smaller rank
	793	CClientClientDHTSizet dhtEdgeHash (edgeHash2Rank, comm);
	794	dhtEdgeHash.computeIndexInfoMapping(edgeHashList);
	795	CClientClientDHTSizet::Index2VectorInfoTypeMap& edgeHash2Info = dhtEdgeHash.getInfoIndexMap();
[900]	796
[1158]	797	// edgeHash2Info = <edgeHash, < rank1, idxGlo, rank2, idxGlo>>
	798	int edgeCount = 0;
[924]	799	for (CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator it = edgeHash2Info.begin(); it != edgeHash2Info.end(); ++it)
	800	{
	801	vector <size_t> edgeInfo = it->second;
[1158]	802	if (edgeInfo[0] == mpiRank)
	803	{
	804	++edgeCount;
	805	}
[924]	806	}
[879]	807
[1507]	808	unsigned long edgeStart, nbEdges;
[1158]	809	MPI_Scan(&edgeCount, &edgeStart, 1, MPI_UNSIGNED_LONG, MPI_SUM, comm);
	810	int nEdges = edgeStart;
	811	MPI_Bcast(&nEdges, 1, MPI_UNSIGNED_LONG, mpiSize-1, comm);
	812	nbEdgesGlo = nEdges;
[924]	813
[1158]	814	// edges to be splitted equally between procs
	815	if ( (nbEdgesGlo % mpiSize) == 0)
	816	{
	817	edge_count = nbEdgesGlo/mpiSize;
	818	edge_start = mpiRank*edge_count;
	819	}
	820	else
	821	{
	822	if (mpiRank == (mpiSize - 1) )
	823	{
	824	edge_count = nbEdgesGlo/mpiSize;
	825	edge_start = mpiRank*(nbEdgesGlo/mpiSize + 1);
	826	}
	827	else
	828	{
	829	edge_count = nbEdgesGlo/mpiSize + 1;
	830	edge_start = mpiRank*edge_count;
	831	}
	832	}
	833	CArray<size_t,1> edgeIdxGloList(edge_count);
	834	for (int i = 0; i < edge_count; ++i)
	835	{
	836	edgeIdxGloList(i) = i + edge_start;
	837	}
[924]	838
[1158]	839	CClientClientDHTSizet dhtEdgeIdxGlo2Face (edgeIdxGlo2Face, comm);
[924]	840	CClientClientDHTSizet dhtEdge2Face (edgeIdxGlo2Face, comm);
[1158]	841	dhtEdgeIdxGlo2Face.computeIndexInfoMapping(edgeIdxGloList);
	842	CClientClientDHTSizet::Index2VectorInfoTypeMap& edgeIdxGlo2FaceIdx = dhtEdgeIdxGlo2Face.getInfoIndexMap();
	843	dhtEdge2Face.computeIndexInfoMapping(edgeIdxList);
	844	CClientClientDHTSizet::Index2VectorInfoTypeMap& edgeIdx2FaceIdx = dhtEdge2Face.getInfoIndexMap();
[924]	845
[1158]	846
	847	edge_faces.resize(2, edge_count);
	848	for (int i = 0; i < edge_count; ++i)
	849	{
	850	CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator it = edgeIdxGlo2FaceIdx.find(i + edge_start);
	851	int indexGlo = it->first;
	852	vector<size_t> faces = it->second;
	853	int face1 = faces[0];
	854	edge_faces(0, indexGlo - edge_start) = face1;
	855	if (faces.size() == 2)
	856	{
	857	int face2 = faces[1];
	858	edge_faces(1, indexGlo - edge_start) = face2;
	859	}
	860	else
	861	{
	862	edge_faces(1, indexGlo - edge_start) = -999;
	863	}
	864	}
	865
	866	size_t tmp;
	867	vector<size_t> tmpVec;
	868	for (CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator it = edgeIdx2FaceIdx.begin(); it != edgeIdx2FaceIdx.end(); it++)
	869	{
	870	tmp = it->first;
	871	tmpVec = it->second;
	872	tmp++;
	873	}
	874
	875	CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itFace1, itFace2, itIndex;
	876	face_faces.resize(nvertex, nbFaces_);
[929]	877	for (int nf = 0; nf < nbFaces_; ++nf)
[924]	878	{
	879	for (int nv1 = 0; nv1 < nvertex; ++nv1)
	880	{
	881	int nh1 = 0;
	882	int nv2 = (nv1 < nvertex -1 ) ? (nv1 + 1) : (nv1 + 1 - nvertex); // cyclic rotation
	883	it1 = nodeHash2IdxGlo.find(nodeHashList((nfnvertex + nv1)4 + nh1));
	884	while (it1 == nodeHash2IdxGlo.end())
	885	{
	886	++nh1;
	887	it1 = nodeHash2IdxGlo.find(nodeHashList((nfnvertex + nv1)4 + nh1));
	888	}
	889	int nh2 = 0;
	890	it2 = nodeHash2IdxGlo.find(nodeHashList((nfnvertex + nv2)4 + nh2));
	891	while (it2 == nodeHash2IdxGlo.end())
	892	{
	893	++nh2;
	894	it2 = nodeHash2IdxGlo.find(nodeHashList((nfnvertex + nv1)4 + nh2));
	895	}
	896
[929]	897	if (it1->second[0] != it2->second[0])
[924]	898	{
[929]	899	size_t faceIdxGlo = nbFacesAccum + nf;
	900	size_t edgeHash = hashPairOrdered(it1->second[0], it2->second[0]);
[1158]	901	itEdgeHash = edgeHash2Info.find(edgeHash);
	902	int edgeIdxGlo = (itEdgeHash->second)[1];
	903
	904	if ( (itEdgeHash->second)[0] == mpiRank)
[924]	905	{
[1158]	906	itFace1 = edgeIdx2FaceIdx.find(edgeIdxGlo);
[929]	907	int face1 = itFace1->second[0];
	908	if (itFace1->second.size() == 1)
	909	{
	910	face_faces(nv1, nf) = 999999;
	911	}
	912	else
	913	{
	914	int face2 = itFace1->second[1];
	915	face_faces(nv1, nf) = (faceIdxGlo == face1 ? face2 : face1);
	916	}
	917	} // edge owner
[924]	918	else
	919	{
[1158]	920	itFace1 = edgeIdx2FaceIdx.find(edgeIdxGlo);
[929]	921	int face1 = itFace1->second[0];
[924]	922	int face2 = itFace1->second[1];
	923	face_faces(nv1, nf) = (faceIdxGlo == face1 ? face2 : face1);
[929]	924	} // not an edge owner
	925	} // node1 != node2
[924]	926	else
	927	{
[929]	928	face_faces(nv1, nf) = 999999;
[924]	929	}
	930	}
	931	}
	932	} // edgesAreWritten
	933
	934	// Case (2): nodes have been previously generated
	935	else if (nodesAreWritten)
[879]	936	{
[924]	937	// (2.1) Generating nodeHashList
	938	vector<size_t> hashValues(4);
[929]	939	CArray<size_t,1> nodeHashList(nbFaces_nvertex4);
	940	for (int nf = 0; nf < nbFaces_; ++nf)
[879]	941	{
[924]	942	for (int nv = 0; nv < nvertex; ++nv)
	943	{
	944	hashValues = CMesh::createHashes(bounds_lon(nv, nf), bounds_lat(nv, nf));
	945	for (int nh = 0; nh < 4; ++nh)
	946	nodeHashList((nfnvertex + nv)4 + nh) = hashValues[nh];
	947	}
[879]	948	}
	949
[924]	950	// (2.2) Recuperating node global indexing and saving face_nodes
	951	// Generating edgeHash2Info = <hash, <idx, rank>> and edgeHashList
	952	pNodeGlobalIndex->computeIndexInfoMapping(nodeHashList);
	953	CClientClientDHTSizet::Index2VectorInfoTypeMap& nodeHash2IdxGlo = pNodeGlobalIndex->getInfoIndexMap();
	954	CClientClientDHTSizet::Index2VectorInfoTypeMap edgeHash2Idx;
	955	CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator it1, it2;
[929]	956	CArray<size_t,1> edgeHashList(nbFaces_*nvertex);
[1158]	957	int nEdgeHash = 0;
[929]	958	for (int nf = 0; nf < nbFaces_; ++nf)
[924]	959	{
	960	for (int nv1 = 0; nv1 < nvertex; ++nv1)
	961	{
	962	int nh1 = 0;
	963	int nv2 = (nv1 < nvertex -1 ) ? (nv1 + 1) : (nv1 + 1 - nvertex); // cyclic rotation
	964	it1 = nodeHash2IdxGlo.find(nodeHashList((nfnvertex + nv1)4 + nh1));
	965	while (it1 == nodeHash2IdxGlo.end())
	966	{
	967	++nh1;
	968	it1 = nodeHash2IdxGlo.find(nodeHashList((nfnvertex + nv1)4 + nh1));
	969	}
	970	int nh2 = 0;
	971	it2 = nodeHash2IdxGlo.find(nodeHashList((nfnvertex + nv2)4 + nh2));
	972	while (it2 == nodeHash2IdxGlo.end())
	973	{
	974	++nh2;
	975	it2 = nodeHash2IdxGlo.find(nodeHashList((nfnvertex + nv1)4 + nh2));
	976	}
	977	face_nodes(nv1,nf) = it1->second[0];
	978	size_t edgeHash = hashPairOrdered(it1->second[0], it2->second[0]);
[1158]	979	if (edgeHash2Idx.count(edgeHash) == 0)
	980	{
	981	edgeHash2Idx[edgeHash].push_back(edgeHash);
	982	edgeHash2Idx[edgeHash].push_back(mpiRank);
	983	edgeHashList(nEdgeHash) = edgeHash;
	984	++nEdgeHash;
	985	}
[924]	986	}
	987	}
[1158]	988	edgeHashList.resizeAndPreserve(nEdgeHash);
[879]	989
[1158]	990	// (2.3) Generating global edge indexes
	991	// The ownership criterion: priority of the process with smaller rank
[924]	992	// Maps generated in this step are:
[1158]	993	// edgeIdx2Idx = = <idx, <rankOwner, idx>>
	994	// edgeIdx2IdxGlo = <idxMin, <rankOwner, idxGlo>>
[924]	995
	996	CClientClientDHTSizet dhtEdgeHash(edgeHash2Idx, comm);
	997	dhtEdgeHash.computeIndexInfoMapping(edgeHashList);
	998	CClientClientDHTSizet::Index2VectorInfoTypeMap& edgeHash2Info = dhtEdgeHash.getInfoIndexMap();
[1158]	999	// edgeHash2Info = <hash, [[idx1, rank1], [idx2, rank2], [idx3, rank3]..]>
[924]	1000
[1158]	1001	CClientClientDHTSizet::Index2VectorInfoTypeMap edgeIdx2Idx;
[924]	1002
	1003	for (CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator it = edgeHash2Info.begin(); it != edgeHash2Info.end(); ++it)
	1004	{
[1158]	1005	size_t rankMin = (it->second)[1];
[924]	1006	size_t idx = (it->second)[0];
[1158]	1007
[924]	1008	for (int i = 2; i < (it->second).size();)
	1009	{
[1158]	1010	if ((it->second)[i+1] < rankMin)
[924]	1011	{
[1158]	1012	rankMin = (it->second)[i+1];
[924]	1013	idx = (it->second)[i];
[1158]	1014	(it->second)[i+1] = (it->second)[i-1];
[924]	1015	}
[1158]	1016	i += 2;
	1017	}
	1018	if (edgeIdx2Idx.count(idx) == 0)
	1019	{
	1020	if (mpiRank == rankMin)
[924]	1021	{
[1158]	1022	edgeIdx2Idx[idx].push_back(rankMin);
	1023	edgeIdx2Idx[idx].push_back(idx);
[924]	1024	}
	1025	}
[1158]	1026	}
	1027
[1507]	1028	unsigned long edgeCount = edgeIdx2Idx.size();
	1029	unsigned long edgeStart, nbEdges;
[1158]	1030	MPI_Scan(&edgeCount, &edgeStart, 1, MPI_UNSIGNED_LONG, MPI_SUM, comm);
	1031	int nEdges = edgeStart;
	1032	MPI_Bcast(&nEdges, 1, MPI_UNSIGNED_LONG, mpiSize-1, comm);
	1033	nbEdgesGlo = nEdges;
	1034
	1035	edgeStart -= edgeCount;
	1036	edge_start = edgeStart;
	1037	edge_count = edgeCount;
	1038	CClientClientDHTSizet::Index2VectorInfoTypeMap dummyEdgeMap;
	1039	int count = 0;
	1040
	1041	for (int nf = 0; nf < nbFaces_; ++nf)
	1042	{
	1043	for (int nv1 = 0; nv1 < nvertex; ++nv1)
[924]	1044	{
[1158]	1045	// Getting global indexes of edge's nodes
	1046	int nh1 = 0;
	1047	int nv2 = (nv1 < nvertex -1 ) ? (nv1 + 1) : (nv1 + 1 - nvertex); // cyclic rotation
	1048	it1 = nodeHash2IdxGlo.find(nodeHashList((nfnvertex + nv1)4 + nh1));
	1049	while (it1 == nodeHash2IdxGlo.end())
	1050	{
	1051	++nh1;
	1052	it1 = nodeHash2IdxGlo.find(nodeHashList((nfnvertex + nv1)4 + nh1));
	1053	}
	1054	int nh2 = 0;
	1055	it2 = nodeHash2IdxGlo.find(nodeHashList((nfnvertex + nv2)4 + nh2));
	1056	while (it2 == nodeHash2IdxGlo.end())
	1057	{
	1058	++nh2;
	1059	it2 = nodeHash2IdxGlo.find(nodeHashList((nfnvertex + nv1)4 + nh2));
	1060	}
	1061	size_t nodeIdxGlo1 = it1->second[0];
	1062	size_t nodeIdxGlo2 = it2->second[0];
	1063
	1064	if (nodeIdxGlo1 != nodeIdxGlo2)
	1065	{
	1066	size_t edgeIdx = hashPairOrdered(nodeIdxGlo1, nodeIdxGlo2);
	1067	CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator it = edgeIdx2Idx.find(edgeIdx);
	1068	if (it != edgeIdx2Idx.end())
	1069	{
	1070	if (dummyEdgeMap.count(edgeIdx) == 0)
	1071	{
	1072	dummyEdgeMap[edgeIdx].push_back(edgeIdx);
	1073	(it->second)[1] = edge_start + count;
	1074	++count;
	1075	}
	1076	}
	1077	}
[924]	1078	}
	1079	}
	1080
[1158]	1081	CClientClientDHTSizet dhtEdgeIdx(edgeIdx2Idx, comm);
	1082	dhtEdgeIdx.computeIndexInfoMapping(edgeHashList);
	1083	CClientClientDHTSizet::Index2VectorInfoTypeMap& edgeIdx2IdxGlo = dhtEdgeIdx.getInfoIndexMap();
	1084
	1085	// (2.4) Saving variables: edge_lon, edge_lat, face_edges
[924]	1086	edge_lon.resize(edge_count);
	1087	edge_lat.resize(edge_count);
	1088	edge_nodes.resize(2, edge_count);
[929]	1089	face_edges.resize(nvertex, nbFaces_);
[924]	1090
	1091	CClientClientDHTSizet::Index2VectorInfoTypeMap edgeIdxGlo2Face;
[929]	1092	CArray<size_t,1> edgeIdxGloList(nbFaces_*nvertex);
[924]	1093	size_t iIdx = 0;
	1094
[929]	1095	for (int nf = 0; nf < nbFaces_; ++nf)
[924]	1096	{
	1097	for (int nv1 = 0; nv1 < nvertex; ++nv1)
	1098	{
	1099	// Getting global indexes of edge's nodes
	1100	int nh1 = 0;
	1101	int nv2 = (nv1 < nvertex -1 ) ? (nv1 + 1) : (nv1 + 1 - nvertex); // cyclic rotation
	1102	it1 = nodeHash2IdxGlo.find(nodeHashList((nfnvertex + nv1)4 + nh1));
	1103	while (it1 == nodeHash2IdxGlo.end())
	1104	{
	1105	++nh1;
	1106	it1 = nodeHash2IdxGlo.find(nodeHashList((nfnvertex + nv1)4 + nh1));
	1107	}
	1108	int nh2 = 0;
	1109	it2 = nodeHash2IdxGlo.find(nodeHashList((nfnvertex + nv2)4 + nh2));
	1110	while (it2 == nodeHash2IdxGlo.end())
	1111	{
	1112	++nh2;
	1113	it2 = nodeHash2IdxGlo.find(nodeHashList((nfnvertex + nv1)4 + nh2));
	1114	}
	1115	// Getting edge global index
	1116	size_t nodeIdxGlo1 = it1->second[0];
	1117	size_t nodeIdxGlo2 = it2->second[0];
[1158]	1118	size_t myIdx = hashPairOrdered(nodeIdxGlo1, nodeIdxGlo2);
[929]	1119	if (nodeIdxGlo1 != nodeIdxGlo2)
	1120	{
[1158]	1121	CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itIdx = edgeIdx2IdxGlo.find(myIdx);
	1122	int edgeIdxGlo = (itIdx->second)[1];
[929]	1123	size_t faceIdxGlo = nbFacesAccum + nf;
[924]	1124
[1158]	1125	if (mpiRank == (itIdx->second)[0])
[929]	1126	{
	1127	double edgeLon;
	1128	double diffLon = abs(bounds_lon(nv1, nf) - bounds_lon(nv2, nf));
	1129	if (diffLon < (180.- prec))
	1130	edgeLon = ( bounds_lon(nv1, nf) + bounds_lon(nv2, nf)) * 0.5;
	1131	else if (diffLon > (180.+ prec))
	1132	edgeLon = (bounds_lon(nv1, nf) + bounds_lon(nv2, nf)) * 0.5 -180.;
	1133	else
	1134	edgeLon = 0.;
	1135	edge_lon(edgeIdxGlo - edge_start) = edgeLon;
	1136	edge_lat(edgeIdxGlo - edge_start) = ( bounds_lat(nv1, nf) + bounds_lat(nv2, nf) ) * 0.5;
	1137	edge_nodes(0, edgeIdxGlo - edge_start) = nodeIdxGlo1;
	1138	edge_nodes(1, edgeIdxGlo - edge_start) = nodeIdxGlo2;
	1139	}
	1140	face_edges(nv1,nf) = edgeIdxGlo;
	1141	if (edgeIdxGlo2Face.count(edgeIdxGlo) == 0)
	1142	{
	1143	edgeIdxGloList(iIdx) = edgeIdxGlo;
	1144	++iIdx;
	1145	}
	1146	edgeIdxGlo2Face[edgeIdxGlo].push_back(faceIdxGlo);
	1147	} // nodeIdxGlo1 != nodeIdxGlo2
[924]	1148	else
	1149	{
[929]	1150	face_edges(nv1,nf) = 999999;
[924]	1151	}
	1152	}
	1153	}
	1154	edgeIdxGloList.resizeAndPreserve(iIdx);
	1155
[1158]	1156	// (2.5) Saving remaining variables edge_faces and face_faces
[924]	1157	edge_faces.resize(2, edge_count);
[929]	1158	face_faces.resize(nvertex, nbFaces_);
[924]	1159
	1160	CClientClientDHTSizet dhtEdge2Face (edgeIdxGlo2Face, comm);
	1161	dhtEdge2Face.computeIndexInfoMapping(edgeIdxGloList);
	1162	CClientClientDHTSizet::Index2VectorInfoTypeMap& edgeIdxGlo2FaceIdx = dhtEdge2Face.getInfoIndexMap();
	1163	CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itNodeIdxGlo1, itNodeIdxGlo2;
[1158]	1164	CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itIdx;
[924]	1165
[929]	1166	for (int nf = 0; nf < nbFaces_; ++nf)
[924]	1167	{
	1168	for (int nv1 = 0; nv1 < nvertex; ++nv1)
	1169	{
	1170	// Getting global indexes of edge's nodes
	1171	int nh1 = 0;
	1172	int nv2 = (nv1 < nvertex -1 ) ? (nv1 + 1) : (nv1 + 1 - nvertex); // cyclic rotation
	1173	it1 = nodeHash2IdxGlo.find(nodeHashList((nfnvertex + nv1)4 + nh1));
	1174	while (it1 == nodeHash2IdxGlo.end())
	1175	{
	1176	++nh1;
	1177	it1 = nodeHash2IdxGlo.find(nodeHashList((nfnvertex + nv1)4 + nh1));
	1178	}
	1179	int nh2 = 0;
	1180	it2 = nodeHash2IdxGlo.find(nodeHashList((nfnvertex + nv2)4 + nh2));
	1181	while (it2 == nodeHash2IdxGlo.end())
	1182	{
	1183	++nh2;
	1184	it2 = nodeHash2IdxGlo.find(nodeHashList((nfnvertex + nv1)4 + nh2));
	1185	}
	1186	size_t nodeIdxGlo1 = it1->second[0];
	1187	size_t nodeIdxGlo2 = it2->second[0];
	1188
[1158]	1189	size_t myIdx = hashPairOrdered(nodeIdxGlo1, nodeIdxGlo2);
	1190	itIdx = edgeIdx2IdxGlo.find(myIdx);
[929]	1191	size_t faceIdxGlo = nbFacesAccum + nf;
[1158]	1192	int edgeIdxGlo = (itIdx->second)[1];
[924]	1193
[1158]	1194	if (mpiRank == (itIdx->second)[0])
[924]	1195	{
	1196	it1 = edgeIdxGlo2FaceIdx.find(edgeIdxGlo);
	1197	int face1 = it1->second[0];
	1198	if (it1->second.size() == 1)
	1199	{
	1200	edge_faces(0, edgeIdxGlo - edge_start) = face1;
	1201	edge_faces(1, edgeIdxGlo - edge_start) = -999;
[929]	1202	face_faces(nv1, nf) = 999999;
[924]	1203	}
	1204	else
	1205	{
	1206	int face2 = it1->second[1];
	1207	edge_faces(0, edgeIdxGlo - edge_start) = face1;
	1208	edge_faces(1, edgeIdxGlo - edge_start) = face2;
	1209	face_faces(nv1, nf) = (faceIdxGlo == face1 ? face2 : face1);
	1210	}
	1211	}
	1212	else
	1213	{
	1214	it1 = edgeIdxGlo2FaceIdx.find(edgeIdxGlo);
	1215	int face1 = it1->second[0];
	1216	int face2 = it1->second[1];
	1217	face_faces(nv1, nf) = (faceIdxGlo == face1 ? face2 : face1);
	1218	}
	1219	}
	1220	}
	1221	} // nodesAreWritten
	1222
	1223	// Case (3): Neither nodes nor edges have been previously generated
	1224	else
[879]	1225	{
[924]	1226	// (3.1) Creating a list of hashes for each node and a map nodeHash2Idx <hash, <idx,rank> >
	1227	vector<size_t> hashValues(4);
	1228	CClientClientDHTSizet::Index2VectorInfoTypeMap nodeHash2Idx;
[929]	1229	CArray<size_t,1> nodeHashList(nbFaces_nvertex4);
[924]	1230	size_t iHash = 0;
[929]	1231	for (int nf = 0; nf < nbFaces_; ++nf)
[879]	1232	{
[924]	1233	for (int nv = 0; nv < nvertex; ++nv)
[879]	1234	{
[924]	1235	hashValues = CMesh::createHashes(bounds_lon(nv, nf), bounds_lat(nv, nf));
[1158]	1236	// size_t nodeIndex = generateNodeIndex(hashValues, mpiRank);
	1237	size_t nodeIndex = generateNodeIndex(hashValues);
[924]	1238	for (int nh = 0; nh < 4; ++nh)
	1239	{
	1240	if (nodeHash2Idx.count(hashValues[nh])==0)
	1241	{
	1242	nodeHash2Idx[hashValues[nh]].push_back(nodeIndex);
	1243	nodeHash2Idx[hashValues[nh]].push_back(mpiRank);
	1244	nodeHashList(iHash) = hashValues[nh];
	1245	++iHash;
	1246	}
	1247	}
[879]	1248	}
[924]	1249	}
	1250	nodeHashList.resizeAndPreserve(iHash);
	1251
	1252	// (3.2) Generating global node indexes
[1158]	1253	// The ownership criterion: priority of the process with smaller rank.
	1254	// With any other criterion it is not possible to have consistent node indexing for different number of procs.
[924]	1255	// Maps generated in this step are:
[1158]	1256	// nodeHash2Info = <hash, [[idx, rankMin], [idx, rank1], [idx, rank3]..]>
	1257	// nodeIdx2Idx = <idx, <rankOwner, idx>>
[924]	1258
	1259	CClientClientDHTSizet dhtNodeHash(nodeHash2Idx, comm);
	1260	dhtNodeHash.computeIndexInfoMapping(nodeHashList);
	1261	CClientClientDHTSizet::Index2VectorInfoTypeMap& nodeHash2Info = dhtNodeHash.getInfoIndexMap();
	1262
[1158]	1263	CClientClientDHTSizet::Index2VectorInfoTypeMap nodeIdx2Idx;
	1264	CArray<size_t,1> nodeIdxList(nbFaces_nvertex4);
	1265	size_t nIdx = 0;
[924]	1266
	1267	for (CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator it = nodeHash2Info.begin(); it != nodeHash2Info.end(); ++it)
	1268	{
[1158]	1269	size_t rankMin = (it->second)[1];
[924]	1270	size_t idx = (it->second)[0];
	1271	for (int i = 2; i < (it->second).size();)
[879]	1272	{
[1158]	1273	if ( (it->second)[i+1] < rankMin)
[924]	1274	{
	1275	idx = (it->second)[i];
[1158]	1276	rankMin = (it->second)[i+1];
	1277	(it->second)[i+1] = (it->second)[i-1];
[924]	1278	}
[1158]	1279	i += 2;
	1280	}
	1281	if (nodeIdx2Idx.count(idx) == 0)
	1282	{
	1283	if (mpiRank == rankMin)
[924]	1284	{
[1158]	1285	nodeIdx2Idx[idx].push_back(rankMin);
	1286	nodeIdx2Idx[idx].push_back(idx);
[924]	1287	}
[1158]	1288	nodeIdxList(nIdx) = idx;
	1289	++nIdx;
[879]	1290	}
[1158]	1291	}
	1292
	1293	// CDHTAutoIndexing dhtNodeIdxGlo = CDHTAutoIndexing(nodeIdx2Idx, comm);
	1294	// CDHTAutoIndexing will not give consistent node numbering for varying number of procs. =>
	1295	// Solution: global node indexing by hand.
	1296	// Maps modified in this step:
	1297	// nodeIdx2Idx = <idx, idxGlo>
[1507]	1298	unsigned long nodeCount = nodeIdx2Idx.size();
	1299	unsigned long nodeStart, nbNodes;
[1158]	1300	MPI_Scan(&nodeCount, &nodeStart, 1, MPI_UNSIGNED_LONG, MPI_SUM, comm);
	1301	int nNodes = nodeStart;
	1302	MPI_Bcast(&nNodes, 1, MPI_UNSIGNED_LONG, mpiSize-1, comm);
	1303	nbNodesGlo = nNodes;
	1304
	1305	nodeStart -= nodeCount;
	1306	node_start = nodeStart;
	1307	node_count = nodeCount;
	1308	CClientClientDHTSizet::Index2VectorInfoTypeMap dummyMap; // just a dummy map used to ensure that each node is numbered only once
	1309	size_t count = 0;
	1310
	1311	for (int nf = 0; nf < nbFaces_; ++nf)
	1312	{
	1313	for (int nv = 0; nv < nvertex; ++nv)
[924]	1314	{
[1158]	1315	vector<size_t> hashValues = CMesh::createHashes(bounds_lon(nv, nf), bounds_lat(nv, nf));
	1316	size_t nodeIdx = generateNodeIndex(hashValues);
	1317	CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator it = nodeIdx2Idx.find(nodeIdx);
	1318	if (it != nodeIdx2Idx.end())
	1319	{
	1320	if (dummyMap.count(nodeIdx) == 0)
	1321	{
	1322	dummyMap[nodeIdx].push_back(nodeIdx);
	1323	(it->second)[1] = node_start + count;
	1324	++count;
	1325	}
	1326	}
[924]	1327	}
[879]	1328	}
[1158]	1329	nodeIdxList.resizeAndPreserve(nIdx);
	1330	CClientClientDHTSizet dhtNodeIdx(nodeIdx2Idx, comm);
	1331	dhtNodeIdx.computeIndexInfoMapping(nodeIdxList);
	1332	CClientClientDHTSizet::Index2VectorInfoTypeMap& nodeIdx2IdxGlo = dhtNodeIdx.getInfoIndexMap();
[879]	1333
[924]	1334	// (3.3) Saving node data: node_lon, node_lat, and face_nodes
	1335	// Generating edgeHash2Info = <hash, <idx, rank>> and edgeHashList
[1158]	1336	// nbNodesGlo = dhtNodeIdxGlo.getNbIndexesGlobal();
	1337	// node_count = dhtNodeIdxGlo.getIndexCount();
	1338	// node_start = dhtNodeIdxGlo.getIndexStart();
[924]	1339	node_lon.resize(node_count);
	1340	node_lat.resize(node_count);
	1341	size_t nodeIdxGlo1 = 0;
	1342	size_t nodeIdxGlo2 = 0;
	1343	CClientClientDHTSizet::Index2VectorInfoTypeMap edgeHash2Idx;
	1344	CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itNodeIdxGlo1, itNodeIdxGlo2;
[929]	1345	CArray<size_t,1> edgeHashList(nbFaces_*nvertex);
	1346	size_t nEdgeHash = 0;
[879]	1347
[929]	1348	for (int nf = 0; nf < nbFaces_; ++nf)
[924]	1349	{
	1350	for (int nv1 = 0; nv1 < nvertex; ++nv1)
	1351	{
	1352	int nv2 = (nv1 < nvertex -1 ) ? (nv1 + 1) : (nv1 + 1 - nvertex); // cyclic rotation
	1353	vector<size_t> hashValues1 = CMesh::createHashes(bounds_lon(nv1, nf), bounds_lat(nv1, nf));
	1354	vector<size_t> hashValues2 = CMesh::createHashes(bounds_lon(nv2, nf), bounds_lat(nv2, nf));
[1158]	1355	size_t nodeIdx1 = generateNodeIndex(hashValues1);
	1356	size_t nodeIdx2 = generateNodeIndex(hashValues2);
	1357	CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itNodeIdx1 = nodeIdx2IdxGlo.find(nodeIdx1);
	1358	CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itNodeIdx2 = nodeIdx2IdxGlo.find(nodeIdx2);
	1359	size_t ownerRank = (itNodeIdx1->second)[0];
	1360	nodeIdxGlo1 = (itNodeIdx1->second)[1];
	1361	nodeIdxGlo2 = (itNodeIdx2->second)[1];
[900]	1362
[1158]	1363	if (mpiRank == ownerRank)
[924]	1364	{
	1365	node_lon(nodeIdxGlo1 - node_start) = bounds_lon(nv1, nf);
	1366	node_lat(nodeIdxGlo1 - node_start) = bounds_lat(nv1, nf);
	1367	}
[929]	1368	if (nodeIdxGlo1 != nodeIdxGlo2)
	1369	{
	1370	size_t edgeHash = hashPairOrdered(nodeIdxGlo1, nodeIdxGlo2);
[1158]	1371	edgeHash2Idx[edgeHash].push_back(edgeHash);
[929]	1372	edgeHash2Idx[edgeHash].push_back(mpiRank);
	1373	edgeHashList(nEdgeHash) = edgeHash;
	1374	++nEdgeHash;
	1375	}
[924]	1376	face_nodes(nv1,nf) = nodeIdxGlo1;
	1377	}
	1378	}
[929]	1379	edgeHashList.resizeAndPreserve(nEdgeHash);
[924]	1380
	1381	// (3.4) Generating global edge indexes
	1382	// Maps generated in this step are:
[1158]	1383	// edgeIdx2Idx = = <idx, <rankOwner, idx>>
	1384	// edgeIdx2IdxGlo = <idxMin, <rankOwner, idxGlo>>
[924]	1385
	1386	CClientClientDHTSizet dhtEdgeHash(edgeHash2Idx, comm);
	1387	dhtEdgeHash.computeIndexInfoMapping(edgeHashList);
	1388	CClientClientDHTSizet::Index2VectorInfoTypeMap& edgeHash2Info = dhtEdgeHash.getInfoIndexMap();
	1389	// edgeHash2Info = <hash, [[idx1, rank1], [idx2, rank2], [idx3, rank3]..]>
	1390
[1158]	1391	CClientClientDHTSizet::Index2VectorInfoTypeMap edgeIdx2Idx;
[924]	1392
	1393	for (CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator it = edgeHash2Info.begin(); it != edgeHash2Info.end(); ++it)
	1394	{
[1158]	1395	size_t rankMin = (it->second)[1];
[924]	1396	size_t idx = (it->second)[0];
	1397
	1398	for (int i = 2; i < (it->second).size();)
	1399	{
[1158]	1400	if ((it->second)[i+1] < rankMin)
[924]	1401	{
[1158]	1402	rankMin = (it->second)[i+1];
[924]	1403	idx = (it->second)[i];
[1158]	1404	(it->second)[i+1] = (it->second)[i-1];
[924]	1405	}
[1158]	1406	i += 2;
	1407	}
	1408	if (edgeIdx2Idx.count(idx) == 0)
	1409	{
	1410	if (mpiRank == rankMin)
[924]	1411	{
[1158]	1412	edgeIdx2Idx[idx].push_back(rankMin);
	1413	edgeIdx2Idx[idx].push_back(idx);
[924]	1414	}
	1415	}
[1158]	1416	}
	1417
[1507]	1418	unsigned long edgeCount = edgeIdx2Idx.size();
	1419	unsigned long edgeStart, nbEdges;
[1158]	1420	MPI_Scan(&edgeCount, &edgeStart, 1, MPI_UNSIGNED_LONG, MPI_SUM, comm);
	1421	int nEdges = edgeStart;
	1422	MPI_Bcast(&nEdges, 1, MPI_UNSIGNED_LONG, mpiSize-1, comm);
	1423	nbEdgesGlo = nEdges;
	1424
	1425	edgeStart -= edgeCount;
	1426	edge_start = edgeStart;
	1427	edge_count = edgeCount;
	1428	CClientClientDHTSizet::Index2VectorInfoTypeMap dummyEdgeMap;
	1429	count = 0;
	1430
	1431	for (int nf = 0; nf < nbFaces_; ++nf)
	1432	{
	1433	for (int nv1 = 0; nv1 < nvertex; ++nv1)
[924]	1434	{
[1158]	1435	int nv2 = (nv1 < nvertex -1 ) ? (nv1 + 1) : (nv1 + 1 - nvertex); // cyclic rotation
	1436	vector<size_t> hashValues1 = CMesh::createHashes(bounds_lon(nv1, nf), bounds_lat(nv1, nf));
	1437	vector<size_t> hashValues2 = CMesh::createHashes(bounds_lon(nv2, nf), bounds_lat(nv2, nf));
	1438	size_t nodeIdx1 = generateNodeIndex(hashValues1);
	1439	size_t nodeIdx2 = generateNodeIndex(hashValues2);
	1440	CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itNodeIdx1 = nodeIdx2IdxGlo.find(nodeIdx1);
	1441	CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itNodeIdx2 = nodeIdx2IdxGlo.find(nodeIdx2);
	1442	nodeIdxGlo1 = (itNodeIdx1->second)[1];
	1443	nodeIdxGlo2 = (itNodeIdx2->second)[1];
	1444
	1445	if (nodeIdxGlo1 != nodeIdxGlo2)
	1446	{
	1447	size_t edgeIdx = hashPairOrdered(nodeIdxGlo1, nodeIdxGlo2);
	1448	CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator it = edgeIdx2Idx.find(edgeIdx);
	1449	if (it != edgeIdx2Idx.end())
	1450	{
	1451	if (dummyEdgeMap.count(edgeIdx) == 0)
	1452	{
	1453	dummyEdgeMap[edgeIdx].push_back(edgeIdx);
	1454	(it->second)[1] = edge_start + count;
	1455	++count;
	1456	}
	1457	}
	1458	}
[924]	1459	}
	1460	}
[1158]	1461	CClientClientDHTSizet dhtEdgeIdx(edgeIdx2Idx, comm);
	1462	dhtEdgeIdx.computeIndexInfoMapping(edgeHashList);
	1463	CClientClientDHTSizet::Index2VectorInfoTypeMap& edgeIdx2IdxGlo = dhtEdgeIdx.getInfoIndexMap();
[924]	1464
	1465	// (3.5) Saving variables: edge_lon, edge_lat, face_edges
	1466	// Creating map edgeIdxGlo2Face <idxGlo, face>
[1158]	1467	// nbEdgesGlo = dhtEdgeIdxGlo.getNbIndexesGlobal();
	1468	// edge_count = dhtEdgeIdxGlo.getIndexCount();
	1469	// edge_start = dhtEdgeIdxGlo.getIndexStart();
[924]	1470
	1471	edge_lon.resize(edge_count);
	1472	edge_lat.resize(edge_count);
	1473	edge_nodes.resize(2, edge_count);
[929]	1474	face_edges.resize(nvertex, nbFaces_);
[924]	1475
	1476	CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator it1, it2;
	1477	CClientClientDHTSizet::Index2VectorInfoTypeMap edgeIdxGlo2Face;
[929]	1478	CArray<size_t,1> edgeIdxGloList(nbFaces_*nvertex);
	1479	size_t nEdge = 0;
[924]	1480
[929]	1481	for (int nf = 0; nf < nbFaces_; ++nf)
[924]	1482	{
	1483	for (int nv1 = 0; nv1 < nvertex; ++nv1)
	1484	{
	1485	// Getting global indexes of edge's nodes
	1486	int nv2 = (nv1 < nvertex -1 ) ? (nv1 + 1) : (nv1 + 1 - nvertex); // cyclic rotation
	1487	vector<size_t> hashValues1 = CMesh::createHashes(bounds_lon(nv1, nf), bounds_lat(nv1, nf));
	1488	vector<size_t> hashValues2 = CMesh::createHashes(bounds_lon(nv2, nf), bounds_lat(nv2, nf));
	1489
[1158]	1490	size_t nodeIdx1 = generateNodeIndex(hashValues1);
	1491	size_t nodeIdx2 = generateNodeIndex(hashValues2);
	1492	it1 = nodeIdx2IdxGlo.find(nodeIdx1);
	1493	it2 = nodeIdx2IdxGlo.find(nodeIdx2);
	1494	size_t nodeIdxGlo1 = (it1->second)[1];
	1495	size_t nodeIdxGlo2 = (it2->second)[1];
[924]	1496
[929]	1497	if (nodeIdxGlo1 != nodeIdxGlo2)
	1498	{
[1158]	1499	size_t myIdx = hashPairOrdered(nodeIdxGlo1, nodeIdxGlo2);
	1500	CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itIdx = edgeIdx2IdxGlo.find(myIdx);
	1501	int edgeIdxGlo = (itIdx->second)[1];
[929]	1502	size_t faceIdxGlo = nbFacesAccum + nf;
[924]	1503
[1158]	1504	if (mpiRank == (itIdx->second)[0])
[929]	1505	{
	1506	double edgeLon;
	1507	double diffLon = abs(bounds_lon(nv1, nf) - bounds_lon(nv2, nf));
	1508	if (diffLon < (180.- prec))
	1509	edgeLon = ( bounds_lon(nv1, nf) + bounds_lon(nv2, nf)) * 0.5;
	1510	else if (diffLon > (180.+ prec))
	1511	edgeLon = (bounds_lon(nv1, nf) + bounds_lon(nv2, nf)) * 0.5 -180.;
	1512	else
	1513	edgeLon = 0.;
	1514	edge_lon(edgeIdxGlo - edge_start) = edgeLon;
	1515	edge_lat(edgeIdxGlo-edge_start) = ( bounds_lat(nv1, nf) + bounds_lat(nv2, nf) ) * 0.5;
	1516	edge_nodes(0, edgeIdxGlo - edge_start) = nodeIdxGlo1;
	1517	edge_nodes(1, edgeIdxGlo - edge_start) = nodeIdxGlo2;
	1518	}
	1519	face_edges(nv1,nf) = edgeIdxGlo;
	1520	if (edgeIdxGlo2Face.count(edgeIdxGlo) == 0)
	1521	{
	1522	edgeIdxGloList(nEdge) = edgeIdxGlo;
	1523	++nEdge;
	1524	}
	1525	edgeIdxGlo2Face[edgeIdxGlo].push_back(faceIdxGlo);
	1526	} // nodeIdxGlo1 != nodeIdxGlo2
[924]	1527	else
	1528	{
[929]	1529	face_edges(nv1,nf) = 999999;
[924]	1530	}
	1531	}
	1532	}
[929]	1533	edgeIdxGloList.resizeAndPreserve(nEdge);
[924]	1534
	1535	// (3.6) Saving remaining variables edge_faces and face_faces
	1536	edge_faces.resize(2, edge_count);
[929]	1537	face_faces.resize(nvertex, nbFaces_);
[924]	1538
	1539	CClientClientDHTSizet dhtEdge2Face (edgeIdxGlo2Face, comm);
	1540	dhtEdge2Face.computeIndexInfoMapping(edgeIdxGloList);
	1541	CClientClientDHTSizet::Index2VectorInfoTypeMap& edgeIdxGlo2FaceIdx = dhtEdge2Face.getInfoIndexMap();
	1542
[929]	1543	for (int nf = 0; nf < nbFaces_; ++nf)
[924]	1544	{
	1545	for (int nv1 = 0; nv1 < nvertex; ++nv1)
	1546	{
	1547	// Getting global indexes of edge's nodes
	1548	int nv2 = (nv1 < nvertex -1 ) ? (nv1 + 1) : (nv1 + 1 - nvertex); // cyclic rotation
	1549	vector<size_t> hashValues1 = CMesh::createHashes(bounds_lon(nv1, nf), bounds_lat(nv1, nf));
	1550	vector<size_t> hashValues2 = CMesh::createHashes(bounds_lon(nv2, nf), bounds_lat(nv2, nf));
	1551
[1158]	1552	size_t myNodeIdx1 = generateNodeIndex(hashValues1);
	1553	size_t myNodeIdx2 = generateNodeIndex(hashValues2);
[929]	1554	if (myNodeIdx1 != myNodeIdx2)
	1555	{
[1158]	1556	it1 = nodeIdx2IdxGlo.find(myNodeIdx1);
	1557	it2 = nodeIdx2IdxGlo.find(myNodeIdx2);
	1558	size_t nodeIdxGlo1 = (it1->second)[1];
	1559	size_t nodeIdxGlo2 = (it2->second)[1];
	1560	size_t myIdx = hashPairOrdered(nodeIdxGlo1, nodeIdxGlo2);
	1561	CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itIdx = edgeIdx2IdxGlo.find(myIdx);
	1562	int edgeIdxGlo = (itIdx->second)[1];
[924]	1563
[929]	1564	size_t faceIdxGlo = nbFacesAccum + nf;
[924]	1565
[1158]	1566	if (mpiRank == (itIdx->second)[0])
[924]	1567	{
[929]	1568	it1 = edgeIdxGlo2FaceIdx.find(edgeIdxGlo);
	1569	int face1 = it1->second[0];
	1570	if (it1->second.size() == 1)
	1571	{
	1572	edge_faces(0, edgeIdxGlo - edge_start) = face1;
	1573	edge_faces(1, edgeIdxGlo - edge_start) = -999;
	1574	face_faces(nv1, nf) = 999999;
	1575	}
	1576	else
	1577	{
	1578	size_t face2 = it1->second[1];
	1579	edge_faces(0, edgeIdxGlo - edge_start) = face1;
	1580	edge_faces(1, edgeIdxGlo - edge_start) = face2;
	1581	face_faces(nv1, nf) = (faceIdxGlo == face1 ? face2 : face1);
	1582	}
[1158]	1583	}
[924]	1584	else
	1585	{
[929]	1586	it1 = edgeIdxGlo2FaceIdx.find(edgeIdxGlo);
	1587	int face1 = it1->second[0];
	1588	int face2 = it1->second[1];
[924]	1589	face_faces(nv1, nf) = (faceIdxGlo == face1 ? face2 : face1);
[1158]	1590	}
[929]	1591	} // myNodeIdx1 != myNodeIdx2
[924]	1592	else
[929]	1593	face_faces(nv1, nf) = 999999;
[924]	1594	}
	1595	}
[1158]	1596
[924]	1597	}
	1598	facesAreWritten = true;
[879]	1599	} // nvertex >= 3
	1600
	1601	} // createMeshEpsilon
	1602
[929]	1603	///----------------------------------------------------------------
	1604	/*!
[931]	1605	* \fn void CMesh::getGloNghbFacesNodeType(const MPI_Comm& comm, const CArray<int, 1>& face_idx,
	1606	const CArray<double, 2>& bounds_lon, const CArray<double, 2>& bounds_lat,
	1607	CArray<int, 2>& nghbFaces)
[929]	1608	* Finds neighboring cells of a local domain for node-type of neighbors.
	1609	* \param [in] comm
[931]	1610	* \param [in] face_idx Array with global indexes.
[929]	1611	* \param [in] bounds_lon Array of boundary longitudes.
	1612	* \param [in] bounds_lat Array of boundary latitudes.
	1613	* \param [out] nghbFaces 2D array of storing global indexes of neighboring cells and their owner procs.
	1614	*/
	1615
[931]	1616	void CMesh::getGloNghbFacesNodeType(const MPI_Comm& comm, const CArray<int, 1>& face_idx,
[929]	1617	const CArray<double, 2>& bounds_lon, const CArray<double, 2>& bounds_lat,
	1618	CArray<int, 2>& nghbFaces)
	1619	{
[931]	1620	int nvertex = bounds_lon.rows();
[929]	1621	int nbFaces = bounds_lon.shape()[1];
	1622	nghbFaces.resize(2, nbFaces*10); // some estimate on max number of neighbouring cells
	1623
	1624	int mpiRank, mpiSize;
	1625	MPI_Comm_rank(comm, &mpiRank);
	1626	MPI_Comm_size(comm, &mpiSize);
	1627
	1628	// (1) Generating unique node indexes
	1629	// (1.1) Creating a list of hashes for each node and a map nodeHash2Idx <hash, <idx,rank> >
	1630	vector<size_t> hashValues(4);
	1631	CClientClientDHTSizet::Index2VectorInfoTypeMap nodeHash2Idx;
	1632	CArray<size_t,1> nodeHashList(nbFacesnvertex4);
	1633	size_t iIdx = 0;
	1634	for (int nf = 0; nf < nbFaces; ++nf)
	1635	{
	1636	for (int nv = 0; nv < nvertex; ++nv)
	1637	{
	1638	hashValues = CMesh::createHashes(bounds_lon(nv, nf), bounds_lat(nv, nf));
	1639	size_t nodeIndex = generateNodeIndex(hashValues, mpiRank);
	1640	for (int nh = 0; nh < 4; ++nh)
	1641	{
	1642	if (nodeHash2Idx.count(hashValues[nh])==0)
	1643	{
	1644	nodeHash2Idx[hashValues[nh]].push_back(nodeIndex);
	1645	nodeHash2Idx[hashValues[nh]].push_back(mpiRank);
	1646	nodeHashList(iIdx) = hashValues[nh];
	1647	++iIdx;
	1648	}
	1649	}
	1650	}
	1651	}
	1652	nodeHashList.resizeAndPreserve(iIdx);
	1653
	1654	// (1.2) Generating node indexes
	1655	// The ownership criterion: priority of the process holding the smaller index
	1656	// Maps generated in this step are:
	1657	// nodeHash2Info = <hash, idx1, idx2, idx3....>
	1658	// nodeIdx2IdxMin = <idx, idxMin>
	1659	// idxMin is a unique node identifier
	1660
	1661	CClientClientDHTSizet dhtNodeHash(nodeHash2Idx, comm);
	1662	dhtNodeHash.computeIndexInfoMapping(nodeHashList);
	1663	CClientClientDHTSizet::Index2VectorInfoTypeMap& nodeHash2Info = dhtNodeHash.getInfoIndexMap();
	1664
	1665	CClientClientDHTSizet::Index2VectorInfoTypeMap nodeIdx2IdxMin;
	1666
	1667	for (CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator it = nodeHash2Info.begin(); it != nodeHash2Info.end(); ++it)
	1668	{
	1669	size_t idxMin = (it->second)[0];
	1670	size_t idx = (it->second)[0];
	1671	for (int i = 2; i < (it->second).size();)
	1672	{
	1673	if (mpiRank == (it->second)[i+1])
	1674	{
	1675	idx = (it->second)[i];
	1676	}
	1677	if ((it->second)[i] < idxMin)
	1678	{
	1679	idxMin = (it->second)[i];
	1680	(it->second)[i] = (it->second)[i-2];
	1681	(it->second)[i+1] = (it->second)[i-1];
	1682	}
	1683	i += 2;
	1684	}
	1685	(it->second)[0] = idxMin;
	1686	if (nodeIdx2IdxMin.count(idx) == 0)
	1687	{
	1688	nodeIdx2IdxMin[idx].push_back(idxMin);
	1689	}
	1690	}
	1691
	1692	// (2) Creating maps nodeIdxMin2Face = <nodeIdxMin, [face1, face2, ...]>
	1693	CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator it;
	1694	CClientClientDHTSizet::Index2VectorInfoTypeMap nodeIdxMin2Face;
	1695	CArray<size_t,1> nodeIdxMinList(nbFacesnvertex4);
	1696
	1697	size_t nNode = 0;
	1698
	1699	for (int nf = 0; nf < nbFaces; ++nf)
	1700	{
	1701	for (int nv = 0; nv < nvertex; ++nv)
	1702	{
	1703	vector<size_t> hashValues = CMesh::createHashes(bounds_lon(nv, nf), bounds_lat(nv, nf));
	1704	size_t myNodeIdx = generateNodeIndex(hashValues, mpiRank);
	1705	it = nodeIdx2IdxMin.find(myNodeIdx);
	1706	size_t nodeIdxMin = (it->second)[0];
[931]	1707	size_t faceIdx = face_idx(nf);
[929]	1708	if (nodeIdxMin2Face.count(nodeIdxMin) == 0)
	1709	{
	1710	nodeIdxMinList(nNode) = nodeIdxMin;
	1711	++nNode;
	1712	}
	1713	nodeIdxMin2Face[nodeIdxMin].push_back(faceIdx);
	1714	nodeIdxMin2Face[nodeIdxMin].push_back(mpiRank);
	1715	}
	1716	}
	1717	nodeIdxMinList.resizeAndPreserve(nNode);
	1718
	1719	// (3) Face_face connectivity
	1720
	1721	// nodeIdxMin2Info = <nodeIdxMin, [face1, face2,...]>
	1722	CClientClientDHTSizet dhtNode2Face (nodeIdxMin2Face, comm);
	1723	dhtNode2Face.computeIndexInfoMapping(nodeIdxMinList);
	1724	CClientClientDHTSizet::Index2VectorInfoTypeMap& nodeIdxMin2Info = dhtNode2Face.getInfoIndexMap();
	1725	CClientClientDHTSizet::Index2VectorInfoTypeMap mapFaces; // auxiliar map
	1726
	1727	int nbNghb = 0;
	1728	CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itNode;
	1729
	1730	for (int nf = 0; nf < nbFaces; ++nf)
	1731	{
	1732	for (int nv = 0; nv < nvertex; ++nv)
	1733	{
	1734	vector<size_t> hashValues = CMesh::createHashes(bounds_lon(nv, nf), bounds_lat(nv, nf));
	1735	size_t myNodeIdx = generateNodeIndex(hashValues, mpiRank);
	1736	itNode = nodeIdx2IdxMin.find(myNodeIdx);
	1737	size_t nodeIdxMin = (itNode->second)[0];
	1738
	1739	itNode = nodeIdxMin2Info.find(nodeIdxMin);
	1740	for (int i = 0; i < itNode->second.size();)
	1741	{
	1742	size_t face = itNode->second[i];
	1743	size_t rank = itNode->second[i+1];
	1744	if (rank != mpiRank)
	1745	if (mapFaces.count(face) == 0)
	1746	{
	1747	nghbFaces(0, nbNghb) = face;
	1748	nghbFaces(1, nbNghb) = rank;
	1749	++nbNghb;
	1750	mapFaces[face].push_back(face);
	1751	}
	1752	i += 2;
	1753	}
	1754	}
	1755	}
	1756	nghbFaces.resizeAndPreserve(2, nbNghb);
[931]	1757	} // getGloNghbFacesNodeType
[929]	1758
	1759	///----------------------------------------------------------------
	1760	/*!
[931]	1761	* \fn void CMesh::getGloNghbFacesEdgeType(const MPI_Comm& comm, const CArray<int, 1>& face_idx,
[929]	1762	const CArray<double, 2>& bounds_lon, const CArray<double, 2>& bounds_lat,
	1763	CArray<int, 2>& nghbFaces)
	1764	* Finds neighboring cells of a local domain for edge-type of neighbors.
	1765	* \param [in] comm
[931]	1766	* \param [in] face_idx Array with global indexes.
[929]	1767	* \param [in] bounds_lon Array of boundary longitudes.
	1768	* \param [in] bounds_lat Array of boundary latitudes.
	1769	* \param [out] nghbFaces 2D array of storing global indexes of neighboring cells and their owner procs.
	1770	*/
	1771
[931]	1772	void CMesh::getGloNghbFacesEdgeType(const MPI_Comm& comm, const CArray<int, 1>& face_idx,
[929]	1773	const CArray<double, 2>& bounds_lon, const CArray<double, 2>& bounds_lat,
	1774	CArray<int, 2>& nghbFaces)
	1775	{
[931]	1776	int nvertex = bounds_lon.rows();
[929]	1777	int nbFaces = bounds_lon.shape()[1];
	1778	nghbFaces.resize(2, nbFaces*10); // estimate of max number of neighbouring cells
	1779
	1780	int mpiRank, mpiSize;
	1781	MPI_Comm_rank(comm, &mpiRank);
	1782	MPI_Comm_size(comm, &mpiSize);
	1783
	1784	// (1) Generating unique node indexes
	1785	// (1.1) Creating a list of hashes for each node and a map nodeHash2Idx <hash, <idx,rank> >
	1786	vector<size_t> hashValues(4);
	1787	CClientClientDHTSizet::Index2VectorInfoTypeMap nodeHash2Idx;
	1788	CArray<size_t,1> nodeHashList(nbFacesnvertex4);
	1789	size_t iIdx = 0;
	1790	for (int nf = 0; nf < nbFaces; ++nf)
	1791	{
	1792	for (int nv = 0; nv < nvertex; ++nv)
	1793	{
	1794	hashValues = CMesh::createHashes(bounds_lon(nv, nf), bounds_lat(nv, nf));
	1795	size_t nodeIndex = generateNodeIndex(hashValues, mpiRank);
	1796	for (int nh = 0; nh < 4; ++nh)
	1797	{
	1798	if (nodeHash2Idx.count(hashValues[nh])==0)
	1799	{
	1800	nodeHash2Idx[hashValues[nh]].push_back(nodeIndex);
	1801	nodeHash2Idx[hashValues[nh]].push_back(mpiRank);
	1802	nodeHashList(iIdx) = hashValues[nh];
	1803	++iIdx;
	1804	}
	1805	}
	1806	}
	1807	}
	1808	nodeHashList.resizeAndPreserve(iIdx);
	1809
	1810	// (1.2) Generating node indexes
	1811	// The ownership criterion: priority of the process holding the smaller index
	1812	// Maps generated in this step are:
	1813	// nodeHash2Info = <hash, idx1, idx2, idx3....>
	1814	// nodeIdx2IdxMin = <idx, idxMin>
	1815	// idxMin is a unique node identifier
	1816
	1817	CClientClientDHTSizet dhtNodeHash(nodeHash2Idx, comm);
	1818	dhtNodeHash.computeIndexInfoMapping(nodeHashList);
	1819	CClientClientDHTSizet::Index2VectorInfoTypeMap& nodeHash2Info = dhtNodeHash.getInfoIndexMap();
	1820
	1821	CClientClientDHTSizet::Index2VectorInfoTypeMap nodeIdx2IdxMin;
	1822
	1823	for (CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator it = nodeHash2Info.begin(); it != nodeHash2Info.end(); ++it)
	1824	{
	1825	size_t idxMin = (it->second)[0];
	1826	size_t idx = (it->second)[0];
	1827	for (int i = 2; i < (it->second).size();)
	1828	{
	1829	if (mpiRank == (it->second)[i+1])
	1830	{
	1831	idx = (it->second)[i];
	1832	}
	1833	if ((it->second)[i] < idxMin)
	1834	{
	1835	idxMin = (it->second)[i];
	1836	(it->second)[i] = (it->second)[i-2];
	1837	(it->second)[i+1] = (it->second)[i-1];
	1838	}
	1839	i += 2;
	1840	}
	1841	(it->second)[0] = idxMin;
	1842	if (nodeIdx2IdxMin.count(idx) == 0)
	1843	{
	1844	nodeIdx2IdxMin[idx].push_back(idxMin);
	1845	}
	1846	}
	1847
	1848	// (2) Creating map edgeHash2Face = <edgeHash, [[face1, rank1], [face2, rank2]]>, where rank1 = rank2 = ...
	1849
	1850	CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator it1, it2, it;
	1851	CClientClientDHTSizet::Index2VectorInfoTypeMap edgeHash2Face;
	1852	CArray<size_t,1> edgeHashList(nbFaces*nvertex);
	1853
	1854	size_t nEdge = 0;
	1855
	1856	for (int nf = 0; nf < nbFaces; ++nf)
	1857	{
	1858	for (int nv1 = 0; nv1 < nvertex; ++nv1)
	1859	{
	1860	// Getting indexes of edge's nodes
	1861	int nv2 = (nv1 < nvertex -1 ) ? (nv1 + 1) : (nv1 + 1 - nvertex); // cyclic rotation
	1862	vector<size_t> hashValues1 = CMesh::createHashes(bounds_lon(nv1, nf), bounds_lat(nv1, nf));
	1863	vector<size_t> hashValues2 = CMesh::createHashes(bounds_lon(nv2, nf), bounds_lat(nv2, nf));
	1864	size_t myNodeIdx1 = generateNodeIndex(hashValues1, mpiRank);
	1865	size_t myNodeIdx2 = generateNodeIndex(hashValues2, mpiRank);
	1866	it1 = nodeIdx2IdxMin.find(myNodeIdx1);
	1867	it2 = nodeIdx2IdxMin.find(myNodeIdx2);
	1868	size_t nodeIdxMin1 = (it1->second)[0];
	1869	size_t nodeIdxMin2 = (it2->second)[0];
[931]	1870	size_t faceIdx = face_idx(nf);
[929]	1871
	1872	if (nodeIdxMin1 != nodeIdxMin2)
	1873	{
	1874	size_t edgeHash = hashPairOrdered(nodeIdxMin1, nodeIdxMin2);
	1875	if (edgeHash2Face.count(edgeHash) == 0)
	1876	{
	1877	edgeHashList(nEdge) = edgeHash;
	1878	++nEdge;
	1879	}
	1880	edgeHash2Face[edgeHash].push_back(faceIdx);
	1881	edgeHash2Face[edgeHash].push_back(mpiRank);
	1882	} // nodeIdxMin1 != nodeIdxMin2
	1883	}
	1884	}
	1885	edgeHashList.resizeAndPreserve(nEdge);
	1886
	1887	// (3) Face_face connectivity
	1888
	1889	int nbNghb = 0;
	1890	CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itNode1, itNode2;
	1891
	1892	// edgeHash2Info = <edgeHash, [[face1, rank1], [face2, rank2]]>
	1893	CClientClientDHTSizet dhtEdge2Face (edgeHash2Face, comm);
	1894	dhtEdge2Face.computeIndexInfoMapping(edgeHashList);
	1895	CClientClientDHTSizet::Index2VectorInfoTypeMap& edgeHash2Info = dhtEdge2Face.getInfoIndexMap();
	1896	CClientClientDHTSizet::Index2VectorInfoTypeMap mapFaces; // auxiliar map
	1897
	1898	for (int nf = 0; nf < nbFaces; ++nf)
	1899	{
	1900	for (int nv1 = 0; nv1 < nvertex; ++nv1)
	1901	{
	1902	// Getting indexes of edge's nodes
	1903	int nv2 = (nv1 < nvertex -1 ) ? (nv1 + 1) : (nv1 + 1 - nvertex); // cyclic rotation
	1904	vector<size_t> hashValues1 = CMesh::createHashes(bounds_lon(nv1, nf), bounds_lat(nv1, nf));
	1905	vector<size_t> hashValues2 = CMesh::createHashes(bounds_lon(nv2, nf), bounds_lat(nv2, nf));
	1906
	1907	size_t myNodeIdx1 = generateNodeIndex(hashValues1, mpiRank);
	1908	size_t myNodeIdx2 = generateNodeIndex(hashValues2, mpiRank);
	1909	itNode1 = nodeIdx2IdxMin.find(myNodeIdx1);
	1910	itNode2 = nodeIdx2IdxMin.find(myNodeIdx2);
	1911	size_t nodeIdxMin1 = (itNode1->second)[0];
	1912	size_t nodeIdxMin2 = (itNode2->second)[0];
	1913
	1914	if (nodeIdxMin1 != nodeIdxMin2)
	1915	{
	1916	size_t edgeHash = hashPairOrdered(nodeIdxMin1, nodeIdxMin2);
	1917	it1 = edgeHash2Info.find(edgeHash);
	1918
	1919	for (int i = 0; i < it1->second.size();)
	1920	{
	1921	size_t face = it1->second[i];
	1922	size_t rank = it1->second[i+1];
	1923	if (rank != mpiRank)
	1924	if (mapFaces.count(face) == 0)
	1925	{
	1926	nghbFaces(0, nbNghb) = face;
	1927	nghbFaces(1, nbNghb) = rank;
	1928	++nbNghb;
	1929	mapFaces[face].push_back(face);
	1930	}
	1931	i += 2;
	1932	}
	1933	} // nodeIdxMin1 != nodeIdxMin2
	1934	}
	1935	}
	1936	nghbFaces.resizeAndPreserve(2, nbNghb);
[931]	1937	} // getGloNghbFacesEdgeType
[929]	1938
	1939	///----------------------------------------------------------------
	1940	/*!
[931]	1941	* \fn void getGlobalNghbFaces (const int nghbType, const MPI_Comm& comm, const CArray<int, 1>& face_idx,
	1942	const CArray<double, 2>& bounds_lon, const CArray<double, 2>& bounds_lat,
	1943	CArray<size_t, 1>& nghbFaces)
	1944	* Finds neighboring faces owned by other procs.
[929]	1945	* \param [in] nghbType 0 for faces sharing nodes, otherwise for faces sharing edges.
	1946	* \param [in] comm
[931]	1947	* \param [in] face_idx Array with global indexes.
[929]	1948	* \param [in] bounds_lon Array of boundary longitudes.
	1949	* \param [in] bounds_lat Array of boundary latitudes.
[931]	1950	* \param [out] nghbFaces 2D array containing neighboring faces and owner ranks.
[929]	1951	*/
	1952
[931]	1953	void CMesh::getGlobalNghbFaces(const int nghbType, const MPI_Comm& comm,
	1954	const CArray<int, 1>& face_idx,
	1955	const CArray<double, 2>& bounds_lon, const CArray<double, 2>& bounds_lat,
	1956	CArray<int, 2>& nghbFaces)
[929]	1957	{
	1958	if (nghbType == 0)
[931]	1959	getGloNghbFacesNodeType(comm, face_idx, bounds_lon, bounds_lat, nghbFaces);
[929]	1960	else
[931]	1961	getGloNghbFacesEdgeType(comm, face_idx, bounds_lon, bounds_lat, nghbFaces);
	1962	} // getGlobalNghbFaces
[929]	1963
[931]	1964	///----------------------------------------------------------------
	1965	/*!
	1966	* \fn void getLocalNghbFaces (const int nghbType,
	1967	const CArray<double, 2>& bounds_lon, const CArray<double, 2>& bounds_lat,
	1968	CArray<size_t, 1>& nghbFaces)
	1969	* \param [in] nghbType 0 for faces sharing nodes, otherwise for faces sharing edges.
	1970	* \param [in] bounds_lon Array of boundary longitudes.
	1971	* \param [in] bounds_lat Array of boundary latitudes.
	1972	* \param [out] nghbFaces 1D array containing neighboring faces.
	1973	*/
	1974
	1975	void CMesh::getLocalNghbFaces(const int nghbType, const CArray<int, 1>& face_idx,
	1976	const CArray<double, 2>& bounds_lon, const CArray<double, 2>& bounds_lat,
	1977	CArray<int, 2>& nghbFaces, CArray<int, 1>& nbNghbFaces)
	1978	{
	1979	if (nghbType == 0)
	1980	getLocNghbFacesNodeType(face_idx, bounds_lon, bounds_lat, nghbFaces, nbNghbFaces);
	1981	else
	1982	getLocNghbFacesEdgeType(face_idx, bounds_lon, bounds_lat, nghbFaces, nbNghbFaces);
	1983	} // getLocalNghbFaces
	1984
	1985	///----------------------------------------------------------------
	1986	/*!
	1987	* \fn void getLocNghbFacesNodeType (const CArray<double, 2>& bounds_lon, const CArray<double, 2>& bounds_lat,
	1988	CArray<int, 2>& nghbFaces)
[945]	1989	* \param [in] face_idx Array with local face indexing.
[931]	1990	* \param [in] bounds_lon Array of boundary longitudes.
	1991	* \param [in] bounds_lat Array of boundary latitudes.
	1992	* \param [out] nghbFaces 2D array containing neighboring faces.
	1993	* \param [out] nbNghbFaces Array containing number of neighboring faces.
	1994	*/
	1995
	1996	void CMesh::getLocNghbFacesNodeType (const CArray<int, 1>& face_idx,
	1997	const CArray<double, 2>& bounds_lon, const CArray<double, 2>& bounds_lat,
	1998	CArray<int, 2>& faceToFaces, CArray<int, 1>& nbNghbFaces)
	1999	{
	2000	int nvertex = bounds_lon.rows();
	2001	int nbFaces = bounds_lon.shape()[1];
	2002	int nbNodes = 0;
	2003	nbNghbFaces.resize(nbFaces);
	2004	nbNghbFaces = 0;
	2005
[946]	2006	// nodeToFaces connectivity
	2007	CClientClientDHTSizet::Index2VectorInfoTypeMap nodeToFaces;
[931]	2008	for (int nf = 0; nf < nbFaces; ++nf)
	2009	for (int nv = 0; nv < nvertex; ++nv)
	2010	{
[946]	2011	size_t nodeHash = (CMesh::createHashes(bounds_lon(nv, nf), bounds_lat(nv ,nf)))[0];
	2012	nodeToFaces[nodeHash].push_back(face_idx(nf));
[931]	2013	}
	2014
	2015	// faceToFaces connectivity
[1542]	2016	std::unordered_map <int, int> mapFaces; // mapFaces = < hash(face1, face2), hash> (the mapped value is irrelevant)
[946]	2017	int maxNb = 20; // some assumption on the max possible number of neighboring cells
	2018	faceToFaces.resize(maxNb, nbFaces);
	2019	CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator it;
	2020	for (it = nodeToFaces.begin(); it != nodeToFaces.end(); ++it)
[931]	2021	{
[946]	2022	int size = it->second.size();
	2023	for (int i = 0; i < (size-1); ++i)
[931]	2024	{
[946]	2025	int face1 = it->second[i];
	2026	for (int j = i+1; j < size; ++j)
[931]	2027	{
[946]	2028	int face2 = it->second[j];
	2029	if (face2 != face1)
[931]	2030	{
[946]	2031	int hashFace = hashPairOrdered(face1, face2);
	2032	if (mapFaces.count(hashFace) == 0)
	2033	{
	2034	faceToFaces(nbNghbFaces(face1), face1) = face2;
	2035	faceToFaces(nbNghbFaces(face2), face2) = face1;
	2036	++nbNghbFaces(face1);
	2037	++nbNghbFaces(face2);
	2038	mapFaces[hashFace] = hashFace;
	2039	}
[931]	2040	}
	2041	}
	2042	}
	2043	}
[946]	2044	} //getLocNghbFacesNodeType
[941]	2045
[931]	2046
	2047	///----------------------------------------------------------------
	2048	/*!
	2049	* \fn void getLocNghbFacesEdgeType (const CArray<int, 1>& face_idx,
	2050	* const CArray<double, 2>& bounds_lon, const CArray<double, 2>& bounds_lat,
	2051	* CArray<int, 2>& nghbFaces, CArray<int, 1>& nbNghbFaces)
[945]	2052	* \param [in] face_idx Array with local face indexing.
[931]	2053	* \param [in] bounds_lon Array of boundary longitudes.
	2054	* \param [in] bounds_lat Array of boundary latitudes.
	2055	* \param [out] nghbFaces 2D array containing neighboring faces.
	2056	* \param [out] nbNghbFaces Array containing number of neighboring faces.
	2057	*/
	2058
	2059	void CMesh::getLocNghbFacesEdgeType (const CArray<int, 1>& face_idx,
	2060	const CArray<double, 2>& bounds_lon, const CArray<double, 2>& bounds_lat,
	2061	CArray<int, 2>& faceToFaces, CArray<int, 1>& nbNghbFaces)
	2062	{
	2063	int nvertex = bounds_lon.rows();
	2064	int nbFaces = bounds_lon.shape()[1];
	2065	int nbNodes = 0;
	2066	int nbEdges = 0;
	2067	nbNghbFaces.resize(nbFaces);
	2068	nbNghbFaces = 0;
	2069
	2070	// faceToNodes connectivity
	2071	CArray<double, 2> faceToNodes (nvertex, nbFaces);
	2072
[1542]	2073	std::unordered_map <pairDouble, int, boost::hash<pairDouble> > mapNodes;
[931]	2074
	2075	for (int nf = 0; nf < nbFaces; ++nf)
	2076	for (int nv = 0; nv < nvertex; ++nv)
	2077	{
	2078	if (mapNodes.find(make_pair (bounds_lon(nv, nf), bounds_lat(nv ,nf))) == mapNodes.end())
	2079	{
	2080	mapNodes[make_pair (bounds_lon(nv, nf), bounds_lat(nv, nf))] = nbNodes;
	2081	faceToNodes(nv,nf) = nbNodes ;
	2082	++nbNodes ;
	2083	}
	2084	else
	2085	faceToNodes(nv,nf) = mapNodes[make_pair (bounds_lon(nv, nf), bounds_lat(nv ,nf))];
	2086	}
	2087
	2088	// faceToFaces connectivity
[1542]	2089	std::unordered_map <pairInt, int, boost::hash<pairInt> > mapEdges;
[931]	2090	faceToFaces.resize(nvertex, nbFaces);
	2091	CArray<int, 2> edgeToFaces(2, nbFaces*nvertex); // max possible
	2092
	2093	for (int nf = 0; nf < nbFaces; ++nf)
	2094	{
	2095	for (int nv1 = 0; nv1 < nvertex; ++nv1)
	2096	{
	2097	int nv2 = (nv1 < nvertex -1 ) ? (nv1 + 1) : (nv1 + 1 - nvertex); // cyclic rotation
	2098	int face = face_idx(nf);
	2099	int node1 = faceToNodes(nv1,face);
	2100	int node2 = faceToNodes(nv2,face);
	2101	if (node1 != node2)
	2102	{
	2103	if (mapEdges.find(make_ordered_pair (node1, node2)) == mapEdges.end())
	2104	{
	2105	mapEdges[make_ordered_pair (node1, node2)] = nbEdges;
	2106	edgeToFaces(0,nbEdges) = face;
	2107	++nbEdges;
	2108	}
	2109	else
	2110	{
	2111	int edge = mapEdges[make_ordered_pair (node1, node2)];
	2112	edgeToFaces(1, edge) = face;
	2113	int face1 = face;
	2114	int face2 = edgeToFaces(0,edge);
	2115	faceToFaces(nbNghbFaces(face1), face1) = face2;
	2116	faceToFaces(nbNghbFaces(face2), face2) = face1;
	2117	++nbNghbFaces(face1);
	2118	++nbNghbFaces(face2);
	2119	}
	2120	} // node1 != node2
	2121	} // nv
	2122	} // nf
	2123
	2124	} //getLocNghbFacesEdgeType
	2125
	2126
[879]	2127	} // namespace xios

Note: See TracBrowser for help on using the repository browser.

Download in other formats: