Context Navigation

source: XIOS/dev/branch_yushan_merged/src/node/mesh.cpp @ 1179

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

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

Download in other formats: