- Timestamp:
- 09/09/16 16:57:43 (8 years ago)
- Location:
- XIOS/trunk/src
- Files:
-
- 3 edited
Legend:
- Unmodified
- Added
- Removed
-
XIOS/trunk/src/io/nc4_data_output.cpp
r924 r929 617 617 SuperClassWriter::addAttribute("long_name", StdString("Maps every face to its edges."), &face_edges); 618 618 SuperClassWriter::addAttribute("start_index", 0, &face_edges); 619 SuperClassWriter::addAttribute("_FillValue", 999999, &face_edges); 619 620 dim0.clear(); 620 621 dim0.push_back(dimEdge); … … 633 634 SuperClassWriter::addAttribute("long_name", StdString("Indicates which other faces neighbor each face"), &face_faces); 634 635 SuperClassWriter::addAttribute("start_index", 0, &face_faces); 636 SuperClassWriter::addAttribute("_FillValue", 999999, &face_faces); 635 637 SuperClassWriter::addAttribute("flag_values", -1, &face_faces); 636 638 SuperClassWriter::addAttribute("flag_meanings", StdString("out_of_mesh"), &face_faces); -
XIOS/trunk/src/node/mesh.cpp
r924 r929 14 14 , node_start{0}, node_count{0} 15 15 , edge_start{0}, edge_count{0} 16 , nbFaces {0}, nbNodes{0}, nbEdges{0}16 , nbFaces_{0}, nbNodes_{0}, nbEdges_{0} 17 17 , nodesAreWritten{false}, edgesAreWritten{false}, facesAreWritten{false} 18 18 , node_lon(), node_lat() … … 310 310 if (nvertex == 1) 311 311 { 312 nbNodes = lonvalue.numElements();313 node_lon.resizeAndPreserve(nbNodes );314 node_lat.resizeAndPreserve(nbNodes );315 for (int nn = 0; nn < nbNodes ; ++nn)312 nbNodes_ = lonvalue.numElements(); 313 node_lon.resizeAndPreserve(nbNodes_); 314 node_lat.resizeAndPreserve(nbNodes_); 315 for (int nn = 0; nn < nbNodes_; ++nn) 316 316 { 317 317 if (map_nodes.find(make_pair (lonvalue(nn), latvalue(nn))) == map_nodes.end()) … … 325 325 else if (nvertex == 2) 326 326 { 327 nbEdges = bounds_lon.shape()[1];327 nbEdges_ = bounds_lon.shape()[1]; 328 328 329 329 // Create nodes and edge_node connectivity 330 node_lon.resizeAndPreserve(nbEdges *nvertex); // Max possible number of nodes331 node_lat.resizeAndPreserve(nbEdges *nvertex);332 edge_nodes.resizeAndPreserve(nvertex, nbEdges );333 334 for (int ne = 0; ne < nbEdges ; ++ne)330 node_lon.resizeAndPreserve(nbEdges_*nvertex); // Max possible number of nodes 331 node_lat.resizeAndPreserve(nbEdges_*nvertex); 332 edge_nodes.resizeAndPreserve(nvertex, nbEdges_); 333 334 for (int ne = 0; ne < nbEdges_; ++ne) 335 335 { 336 336 for (int nv = 0; nv < nvertex; ++nv) … … 338 338 if (map_nodes.find(make_pair (bounds_lon(nv, ne), bounds_lat(nv ,ne))) == map_nodes.end()) 339 339 { 340 map_nodes[make_pair (bounds_lon(nv, ne), bounds_lat(nv, ne))] = nbNodes ;341 edge_nodes(nv,ne) = nbNodes ;342 node_lon(nbNodes ) = bounds_lon(nv, ne);343 node_lat(nbNodes ) = bounds_lat(nv, ne);344 ++nbNodes ;340 map_nodes[make_pair (bounds_lon(nv, ne), bounds_lat(nv, ne))] = nbNodes_ ; 341 edge_nodes(nv,ne) = nbNodes_ ; 342 node_lon(nbNodes_) = bounds_lon(nv, ne); 343 node_lat(nbNodes_) = bounds_lat(nv, ne); 344 ++nbNodes_ ; 345 345 } 346 346 else … … 348 348 } 349 349 } 350 node_lon.resizeAndPreserve(nbNodes );351 node_lat.resizeAndPreserve(nbNodes );350 node_lon.resizeAndPreserve(nbNodes_); 351 node_lat.resizeAndPreserve(nbNodes_); 352 352 353 353 // Create edges 354 edge_lon.resizeAndPreserve(nbEdges );355 edge_lat.resizeAndPreserve(nbEdges );356 357 for (int ne = 0; ne < nbEdges ; ++ne)354 edge_lon.resizeAndPreserve(nbEdges_); 355 edge_lat.resizeAndPreserve(nbEdges_); 356 357 for (int ne = 0; ne < nbEdges_; ++ne) 358 358 { 359 359 if (map_edges.find(make_ordered_pair (edge_nodes(0,ne), edge_nodes(1,ne))) == map_edges.end()) … … 369 369 else 370 370 { 371 nbFaces = bounds_lon.shape()[1];371 nbFaces_ = bounds_lon.shape()[1]; 372 372 373 373 // Create nodes and face_node connectivity 374 node_lon.resizeAndPreserve(nbFaces *nvertex); // Max possible number of nodes375 node_lat.resizeAndPreserve(nbFaces *nvertex);376 face_nodes.resize(nvertex, nbFaces );374 node_lon.resizeAndPreserve(nbFaces_*nvertex); // Max possible number of nodes 375 node_lat.resizeAndPreserve(nbFaces_*nvertex); 376 face_nodes.resize(nvertex, nbFaces_); 377 377 378 for (int nf = 0; nf < nbFaces ; ++nf)378 for (int nf = 0; nf < nbFaces_; ++nf) 379 379 { 380 380 for (int nv = 0; nv < nvertex; ++nv) … … 382 382 if (map_nodes.find(make_pair (bounds_lon(nv, nf), bounds_lat(nv ,nf))) == map_nodes.end()) 383 383 { 384 map_nodes[make_pair (bounds_lon(nv, nf), bounds_lat(nv, nf))] = nbNodes ;385 face_nodes(nv,nf) = nbNodes ;386 node_lon(nbNodes ) = bounds_lon(nv, nf);387 node_lat(nbNodes ) = bounds_lat(nv ,nf);388 ++nbNodes ;384 map_nodes[make_pair (bounds_lon(nv, nf), bounds_lat(nv, nf))] = nbNodes_ ; 385 face_nodes(nv,nf) = nbNodes_ ; 386 node_lon(nbNodes_) = bounds_lon(nv, nf); 387 node_lat(nbNodes_) = bounds_lat(nv ,nf); 388 ++nbNodes_ ; 389 389 } 390 390 else … … 394 394 } 395 395 } 396 node_lon.resizeAndPreserve(nbNodes );397 node_lat.resizeAndPreserve(nbNodes );396 node_lon.resizeAndPreserve(nbNodes_); 397 node_lat.resizeAndPreserve(nbNodes_); 398 398 399 399 // Create edges and edge_nodes connectivity 400 edge_lon.resizeAndPreserve(nbFaces *nvertex); // Max possible number of edges401 edge_lat.resizeAndPreserve(nbFaces *nvertex);402 edge_nodes.resizeAndPreserve(2, nbFaces *nvertex);403 edge_faces.resize(2, nbFaces *nvertex);404 face_edges.resize(nvertex, nbFaces );405 face_faces.resize(nvertex, nbFaces );406 407 vector<int> countEdges(nbFaces *nvertex); // needed in case if edges have been already generated408 vector<int> countFaces(nbFaces );409 countEdges.assign(nbFaces *nvertex, 0);410 countFaces.assign(nbFaces , 0);400 edge_lon.resizeAndPreserve(nbFaces_*nvertex); // Max possible number of edges 401 edge_lat.resizeAndPreserve(nbFaces_*nvertex); 402 edge_nodes.resizeAndPreserve(2, nbFaces_*nvertex); 403 edge_faces.resize(2, nbFaces_*nvertex); 404 face_edges.resize(nvertex, nbFaces_); 405 face_faces.resize(nvertex, nbFaces_); 406 407 vector<int> countEdges(nbFaces_*nvertex); // needed in case if edges have been already generated 408 vector<int> countFaces(nbFaces_); 409 countEdges.assign(nbFaces_*nvertex, 0); 410 countFaces.assign(nbFaces_, 0); 411 411 int edge; 412 for (int nf = 0; nf < nbFaces ; ++nf)412 for (int nf = 0; nf < nbFaces_; ++nf) 413 413 { 414 414 for (int nv1 = 0; nv1 < nvertex; ++nv1) … … 418 418 if (map_edges.find(make_ordered_pair (face_nodes(nv1,nf), face_nodes(nv2,nf))) == map_edges.end()) 419 419 { 420 map_edges[make_ordered_pair (face_nodes(nv1,nf), face_nodes(nv2,nf))] = nbEdges ;420 map_edges[make_ordered_pair (face_nodes(nv1,nf), face_nodes(nv2,nf))] = nbEdges_ ; 421 421 face_edges(nv1,nf) = map_edges[make_ordered_pair (face_nodes(nv1,nf), face_nodes(nv2,nf))]; 422 edge_faces(0,nbEdges ) = nf;423 edge_faces(1,nbEdges ) = -999;424 face_faces(nv1,nf) = -1;425 edge_nodes(Range::all(),nbEdges ) = face_nodes(nv1,nf), face_nodes(nv2,nf);426 edge_lon(nbEdges ) = ( abs( node_lon(face_nodes(nv1,nf)) - node_lon(face_nodes(nv2,nf))) < 180.) ?422 edge_faces(0,nbEdges_) = nf; 423 edge_faces(1,nbEdges_) = -999; 424 face_faces(nv1,nf) = 999999; 425 edge_nodes(Range::all(),nbEdges_) = face_nodes(nv1,nf), face_nodes(nv2,nf); 426 edge_lon(nbEdges_) = ( abs( node_lon(face_nodes(nv1,nf)) - node_lon(face_nodes(nv2,nf))) < 180.) ? 427 427 (( node_lon(face_nodes(nv1,nf)) + node_lon(face_nodes(nv2,nf))) * 0.5) : 428 428 (( node_lon(face_nodes(nv1,nf)) + node_lon(face_nodes(nv2,nf))) * 0.5 -180.); 429 edge_lat(nbEdges ) = ( node_lat(face_nodes(nv1,nf)) + node_lat(face_nodes(nv2,nf)) ) * 0.5;430 ++nbEdges ;429 edge_lat(nbEdges_) = ( node_lat(face_nodes(nv1,nf)) + node_lat(face_nodes(nv2,nf)) ) * 0.5; 430 ++nbEdges_; 431 431 } 432 432 else … … 439 439 if (countEdges[edge]==0) 440 440 { 441 face_faces(nv1,nf) = -1;441 face_faces(nv1,nf) = 999999; 442 442 } 443 443 else … … 465 465 } 466 466 } 467 edge_nodes.resizeAndPreserve(2, nbEdges );468 edge_faces.resizeAndPreserve(2, nbEdges );469 edge_lon.resizeAndPreserve(nbEdges );470 edge_lat.resizeAndPreserve(nbEdges );467 edge_nodes.resizeAndPreserve(2, nbEdges_); 468 edge_faces.resizeAndPreserve(2, nbEdges_); 469 edge_lon.resizeAndPreserve(nbEdges_); 470 edge_lat.resizeAndPreserve(nbEdges_); 471 471 472 472 // Create faces 473 face_lon.resize(nbFaces );474 face_lat.resize(nbFaces );473 face_lon.resize(nbFaces_); 474 face_lat.resize(nbFaces_); 475 475 face_lon = lonvalue; 476 476 face_lat = latvalue; … … 502 502 MPI_Comm_rank(comm, &mpiRank); 503 503 MPI_Comm_size(comm, &mpiSize); 504 double prec = 1e-11; // used in calculations of edge_lon/lat 504 505 505 506 if (nvertex == 1) 506 507 { 507 nbNodes = lonvalue.numElements();508 node_lon.resize(nbNodes );509 node_lat.resize(nbNodes );508 nbNodes_ = lonvalue.numElements(); 509 node_lon.resize(nbNodes_); 510 node_lat.resize(nbNodes_); 510 511 node_lon = lonvalue; 511 512 node_lat = latvalue; … … 514 515 vector<size_t> hashValues(4); 515 516 CClientClientDHTSizet::Index2VectorInfoTypeMap nodeHash2IdxGlo; 516 for (size_t nn = 0; nn < nbNodes ; ++nn)517 for (size_t nn = 0; nn < nbNodes_; ++nn) 517 518 { 518 519 hashValues = CMesh::createHashes(lonvalue(nn), latvalue(nn)); 519 520 for (size_t nh = 0; nh < 4; ++nh) 520 521 { 521 nodeHash2IdxGlo[hashValues[nh]].push_back(mpiRank*nbNodes + nn);522 nodeHash2IdxGlo[hashValues[nh]].push_back(mpiRank*nbNodes_ + nn); 522 523 } 523 524 } … … 528 529 else if (nvertex == 2) 529 530 { 530 nbEdges = bounds_lon.shape()[1];531 edge_lon.resize(nbEdges );532 edge_lat.resize(nbEdges );531 nbEdges_ = bounds_lon.shape()[1]; 532 edge_lon.resize(nbEdges_); 533 edge_lat.resize(nbEdges_); 533 534 edge_lon = lonvalue; 534 535 edge_lat = latvalue; 535 edge_nodes.resize(nvertex, nbEdges); 536 edge_nodes.resize(nvertex, nbEdges_); 537 538 // For determining the global edge index 539 size_t nbEdgesOnProc = nbEdges_; 540 size_t nbEdgesAccum; 541 MPI_Scan(&nbEdgesOnProc, &nbEdgesAccum, 1, MPI_UNSIGNED_LONG, MPI_SUM, comm); 542 nbEdgesAccum -= nbEdges_; 543 536 544 CClientClientDHTSizet::Index2VectorInfoTypeMap edgeHash2IdxGlo; 537 545 CClientClientDHTSizet::Index2VectorInfoTypeMap edgeHash2Idx; … … 541 549 { 542 550 vector<size_t> hashValues(4); 543 CArray<size_t,1> nodeHashList(nbEdges *nvertex*4);544 for (int ne = 0; ne < nbEdges ; ++ne) // size_t doesn't work with CArray<double, 2>551 CArray<size_t,1> nodeHashList(nbEdges_*nvertex*4); 552 for (int ne = 0; ne < nbEdges_; ++ne) // size_t doesn't work with CArray<double, 2> 545 553 { 546 554 for (int nv = 0; nv < nvertex; ++nv) … … 559 567 CClientClientDHTSizet::Index2VectorInfoTypeMap& nodeHash2IdxGlo = pNodeGlobalIndex->getInfoIndexMap(); 560 568 CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator it; 561 vector <size_t> edgeNodes;562 for (int ne = 0; ne < nbEdges ; ++ne)569 size_t nodeIdxGlo1, nodeIdxGlo2; 570 for (int ne = 0; ne < nbEdges_; ++ne) 563 571 { 564 572 for (int nv = 0; nv < nvertex; ++nv) … … 574 582 } 575 583 edge_nodes(nv,ne) = it->second[0]; 576 edgeNodes.push_back(it->second[0]); 577 } 578 edgeHash2IdxGlo[ hashPairOrdered(edgeNodes[0], edgeNodes[1]) ].push_back(nbEdges*mpiRank + ne); 584 if (nv ==0) 585 nodeIdxGlo1 = it->second[0]; 586 else 587 nodeIdxGlo2 = it->second[0]; 588 } 589 size_t edgeIdxGlo = nbEdgesAccum + ne; 590 edgeHash2IdxGlo[ hashPairOrdered(nodeIdxGlo1, nodeIdxGlo2) ].push_back(edgeIdxGlo); 579 591 } 580 592 } // nodesAreWritten 593 581 594 582 595 // Case (2): node indexes have not been generated previously … … 586 599 vector<size_t> hashValues(4); 587 600 CClientClientDHTSizet::Index2VectorInfoTypeMap nodeHash2Idx; 588 CArray<size_t,1> nodeHashList(nbEdges *nvertex*4);589 for (int ne = 0; ne < nbEdges ; ++ne)601 CArray<size_t,1> nodeHashList(nbEdges_*nvertex*4); 602 for (int ne = 0; ne < nbEdges_; ++ne) 590 603 { 591 604 for (int nv = 0; nv < nvertex; ++nv) … … 613 626 CClientClientDHTSizet::Index2VectorInfoTypeMap nodeIdx2IdxMin; 614 627 CClientClientDHTSizet::Index2VectorInfoTypeMap nodeIdx2IdxGlo; 615 CArray<size_t,1> nodeIdxMinList(nbEdges *nvertex);628 CArray<size_t,1> nodeIdxMinList(nbEdges_*nvertex); 616 629 617 630 CClientClientDHTSizet dhtNodeHash(nodeHash2Idx, comm); … … 667 680 size_t idxGlo = 0; 668 681 669 for (int ne = 0; ne < nbEdges ; ++ne)682 for (int ne = 0; ne < nbEdges_; ++ne) 670 683 { 671 684 for (int nv = 0; nv < nvertex; ++nv) … … 695 708 edgeNodes.push_back(idxGlo); 696 709 } 697 edgeHash2IdxGlo[ hashPairOrdered(edgeNodes[0], edgeNodes[1]) ].push_back(nbEdges*mpiRank + ne); 710 if (edgeNodes[0] != edgeNodes[1]) 711 { 712 size_t edgeIdxGlo = nbEdgesAccum + ne; 713 edgeHash2IdxGlo[ hashPairOrdered(edgeNodes[0], edgeNodes[1]) ].push_back(edgeIdxGlo); 714 } 698 715 edgeNodes.clear(); 699 716 } … … 707 724 else 708 725 { 709 nbFaces = bounds_lon.shape()[1];710 face_lon.resize(nbFaces );711 face_lat.resize(nbFaces );726 nbFaces_ = bounds_lon.shape()[1]; 727 face_lon.resize(nbFaces_); 728 face_lat.resize(nbFaces_); 712 729 face_lon = lonvalue; 713 730 face_lat = latvalue; 714 face_nodes.resize(nvertex, nbFaces); 715 face_edges.resize(nvertex, nbFaces); 731 face_nodes.resize(nvertex, nbFaces_); 732 face_edges.resize(nvertex, nbFaces_); 733 734 // For determining the global face index 735 size_t nbFacesOnProc = nbFaces_; 736 size_t nbFacesAccum; 737 MPI_Scan(&nbFacesOnProc, &nbFacesAccum, 1, MPI_UNSIGNED_LONG, MPI_SUM, comm); 738 nbFacesAccum -= nbFaces_; 716 739 717 740 // Case (1): edges have been previously generated … … 720 743 // (1.1) Recuperating node global indexing and saving face_nodes 721 744 vector<size_t> hashValues(4); 722 CArray<size_t,1> nodeHashList(nbFaces *nvertex*4);723 for (int nf = 0; nf < nbFaces ; ++nf)745 CArray<size_t,1> nodeHashList(nbFaces_*nvertex*4); 746 for (int nf = 0; nf < nbFaces_; ++nf) 724 747 { 725 748 for (int nv = 0; nv < nvertex; ++nv) … … 733 756 CClientClientDHTSizet::Index2VectorInfoTypeMap& nodeHash2IdxGlo = pNodeGlobalIndex->getInfoIndexMap(); 734 757 CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator it1, it2; 735 CArray<size_t,1> edgeHashList(nbFaces*nvertex); 736 for (int nf = 0; nf < nbFaces; ++nf) 758 CArray<size_t,1> edgeHashList(nbFaces_*nvertex); 759 size_t nEdge = 0; 760 for (int nf = 0; nf < nbFaces_; ++nf) 737 761 { 738 762 for (int nv1 = 0; nv1 < nvertex; ++nv1) … … 754 778 } 755 779 face_nodes(nv1,nf) = it1->second[0]; 756 edgeHashList(nf*nvertex + nv1) = hashPairOrdered(it1->second[0], it2->second[0]); 757 758 } 759 } 780 if (it1->second[0] != it2->second[0]) 781 { 782 edgeHashList(nEdge) = hashPairOrdered(it1->second[0], it2->second[0]); 783 ++nEdge; 784 } 785 } 786 } 787 edgeHashList.resizeAndPreserve(nEdge); 760 788 761 789 // (1.2) Recuperating edge global indexing and saving face_edges … … 764 792 CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itEdgeHash; 765 793 CClientClientDHTSizet::Index2VectorInfoTypeMap edgeHash2Rank; 766 767 794 CClientClientDHTSizet::Index2VectorInfoTypeMap edgeIdxGlo2Face; 768 CArray<size_t,1> edgeIdxGloList(nbFaces *nvertex);795 CArray<size_t,1> edgeIdxGloList(nbFaces_*nvertex); 769 796 size_t iIdx = 0; 770 797 771 for (int nf = 0; nf < nbFaces ; ++nf)798 for (int nf = 0; nf < nbFaces_; ++nf) 772 799 { 773 800 for (int nv1 = 0; nv1 < nvertex; ++nv1) … … 788 815 it2 = nodeHash2IdxGlo.find(nodeHashList((nf*nvertex + nv1)*4 + nh2)); 789 816 } 790 size_t faceIdxGlo = mpiRank * nbFaces + nf; 791 size_t edgeHash = hashPairOrdered(it1->second[0], it2->second[0]); 792 itEdgeHash = edgeHash2IdxGlo.find(edgeHash); 793 size_t edgeIdxGlo = itEdgeHash->second[0]; 794 face_edges(nv1,nf) = edgeIdxGlo; 795 if (edgeIdxGlo2Face.count(edgeIdxGlo) == 0) 796 { 797 edgeIdxGloList(iIdx) = edgeIdxGlo; 798 ++iIdx; 799 } 800 edgeIdxGlo2Face[edgeIdxGlo].push_back(faceIdxGlo); 801 edgeHash2Rank[edgeHash].push_back(itEdgeHash->first); 802 edgeHash2Rank[edgeHash].push_back(mpiRank); 817 if (it1->second[0] != it2->second[0]) 818 { 819 size_t faceIdxGlo = nbFacesAccum + nf; 820 size_t edgeHash = hashPairOrdered(it1->second[0], it2->second[0]); 821 itEdgeHash = edgeHash2IdxGlo.find(edgeHash); 822 size_t edgeIdxGlo = itEdgeHash->second[0]; 823 face_edges(nv1,nf) = edgeIdxGlo; 824 if (edgeIdxGlo2Face.count(edgeIdxGlo) == 0) 825 { 826 edgeIdxGloList(iIdx) = edgeIdxGlo; 827 ++iIdx; 828 } 829 edgeIdxGlo2Face[edgeIdxGlo].push_back(faceIdxGlo); 830 edgeHash2Rank[edgeHash].push_back(itEdgeHash->first); 831 edgeHash2Rank[edgeHash].push_back(mpiRank); 832 } 833 else 834 { 835 face_edges(nv1,nf) = 999999; 836 } 803 837 } 804 838 } … … 812 846 dhtEdgeHash.computeIndexInfoMapping(edgeHashList); 813 847 CClientClientDHTSizet::Index2VectorInfoTypeMap& edgeHash2Info = dhtEdgeHash.getInfoIndexMap(); 814 CClientClientDHTSizet::Index2VectorInfoTypeMap edgeHashMin2IdxGlo; // map is needed purely for counting848 CClientClientDHTSizet::Index2VectorInfoTypeMap edgeHashMin2IdxGlo; 815 849 816 850 // edgeHash2Info = <edgeHash, <idxGlo, rank1, idxGlo, rank2>> 817 // edgeHashMin2IdxGlo edgeHash2Info= <edgeHashMin, idxGlo>851 // edgeHashMin2IdxGlo = <edgeHashMin, idxGlo> 818 852 for (CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator it = edgeHash2Info.begin(); it != edgeHash2Info.end(); ++it) 819 853 { … … 831 865 832 866 edge_faces.resize(2, edge_count); 833 face_faces.resize(nvertex, nbFaces );867 face_faces.resize(nvertex, nbFaces_); 834 868 835 869 CClientClientDHTSizet dhtEdge2Face (edgeIdxGlo2Face, comm); … … 838 872 CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itFace1, itFace2; 839 873 840 for (int nf = 0; nf < nbFaces ; ++nf)874 for (int nf = 0; nf < nbFaces_; ++nf) 841 875 { 842 876 for (int nv1 = 0; nv1 < nvertex; ++nv1) … … 857 891 it2 = nodeHash2IdxGlo.find(nodeHashList((nf*nvertex + nv1)*4 + nh2)); 858 892 } 859 size_t faceIdxGlo = mpiRank * nbFaces + nf; 860 size_t edgeHash = hashPairOrdered(it1->second[0], it2->second[0]); 861 itEdgeHash = edgeHash2Info.find(edgeHash); 862 size_t edgeOwnerRank = itEdgeHash->second[1]; 863 int edgeIdxGlo = itEdgeHash->second[0]; 864 865 if (mpiRank == edgeOwnerRank) 866 { 867 itFace1 = edgeIdxGlo2FaceIdx.find(edgeIdxGlo); 868 int face1 = itFace1->second[0]; 869 if (itFace1->second.size() == 1) 893 894 if (it1->second[0] != it2->second[0]) 895 { 896 size_t faceIdxGlo = nbFacesAccum + nf; 897 size_t edgeHash = hashPairOrdered(it1->second[0], it2->second[0]); 898 itEdgeHash = edgeHash2IdxGlo.find(edgeHash); 899 if (itEdgeHash != edgeHashMin2IdxGlo.end()) 870 900 { 871 edge_faces(0, edgeIdxGlo - edge_start) = face1; 872 edge_faces(1, edgeIdxGlo - edge_start) = -999; 873 face_faces(nv1, nf) = -1; 874 } 901 int edgeIdxGlo = itEdgeHash->second[0]; 902 itFace1 = edgeIdxGlo2FaceIdx.find(edgeIdxGlo); 903 int face1 = itFace1->second[0]; 904 if (itFace1->second.size() == 1) 905 { 906 edge_faces(0, edgeIdxGlo - edge_start) = face1; 907 edge_faces(1, edgeIdxGlo - edge_start) = -999; 908 face_faces(nv1, nf) = 999999; 909 } 910 else 911 { 912 int face2 = itFace1->second[1]; 913 edge_faces(0, edgeIdxGlo - edge_start) = face1; 914 edge_faces(1, edgeIdxGlo - edge_start) = face2; 915 face_faces(nv1, nf) = (faceIdxGlo == face1 ? face2 : face1); 916 } 917 } // edge owner 875 918 else 876 919 { 920 itEdgeHash = edgeHashMin2IdxGlo.find(edgeHash); 921 size_t edgeIdxGlo = itEdgeHash->second[0]; 922 itFace1 = edgeIdxGlo2FaceIdx.find(edgeIdxGlo); 923 int face1 = itFace1->second[0]; 877 924 int face2 = itFace1->second[1]; 878 edge_faces(0, edgeIdxGlo - edge_start) = face1;879 edge_faces(1, edgeIdxGlo - edge_start) = face2;880 925 face_faces(nv1, nf) = (faceIdxGlo == face1 ? face2 : face1); 881 } 882 } 926 } // not an edge owner 927 } // node1 != node2 883 928 else 884 929 { 885 itFace1 = edgeIdxGlo2FaceIdx.find(edgeIdxGlo); 886 int face1 = itFace1->second[0]; 887 int face2 = itFace1->second[1]; 888 face_faces(nv1, nf) = (faceIdxGlo == face1 ? face2 : face1); 930 face_faces(nv1, nf) = 999999; 889 931 } 890 932 } … … 897 939 // (2.1) Generating nodeHashList 898 940 vector<size_t> hashValues(4); 899 CArray<size_t,1> nodeHashList(nbFaces *nvertex*4);900 for (int nf = 0; nf < nbFaces ; ++nf)941 CArray<size_t,1> nodeHashList(nbFaces_*nvertex*4); 942 for (int nf = 0; nf < nbFaces_; ++nf) 901 943 { 902 944 for (int nv = 0; nv < nvertex; ++nv) … … 914 956 CClientClientDHTSizet::Index2VectorInfoTypeMap edgeHash2Idx; 915 957 CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator it1, it2; 916 CArray<size_t,1> edgeHashList(nbFaces *nvertex);917 for (int nf = 0; nf < nbFaces ; ++nf)958 CArray<size_t,1> edgeHashList(nbFaces_*nvertex); 959 for (int nf = 0; nf < nbFaces_; ++nf) 918 960 { 919 961 for (int nv1 = 0; nv1 < nvertex; ++nv1) … … 955 997 CClientClientDHTSizet::Index2VectorInfoTypeMap edgeIdx2IdxMin; 956 998 CClientClientDHTSizet::Index2VectorInfoTypeMap edgeIdx2IdxGlo; 957 CArray<size_t,1> edgeIdxMinList(nbFaces *nvertex);999 CArray<size_t,1> edgeIdxMinList(nbFaces_*nvertex); 958 1000 size_t iIdxMin = 0; 959 1001 … … 1000 1042 edge_lat.resize(edge_count); 1001 1043 edge_nodes.resize(2, edge_count); 1002 face_edges.resize(nvertex, nbFaces );1044 face_edges.resize(nvertex, nbFaces_); 1003 1045 1004 1046 CClientClientDHTSizet::Index2VectorInfoTypeMap edgeIdxGlo2Face; 1005 CArray<size_t,1> edgeIdxGloList(nbFaces *nvertex);1047 CArray<size_t,1> edgeIdxGloList(nbFaces_*nvertex); 1006 1048 size_t iIdx = 0; 1007 1049 1008 for (int nf = 0; nf < nbFaces ; ++nf)1050 for (int nf = 0; nf < nbFaces_; ++nf) 1009 1051 { 1010 1052 for (int nv1 = 0; nv1 < nvertex; ++nv1) … … 1030 1072 size_t nodeIdxGlo2 = it2->second[0]; 1031 1073 size_t myIdx = generateEdgeIndex(nodeIdxGlo1, nodeIdxGlo2, mpiRank); 1032 CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itIdx = edgeIdx2IdxMin.find(myIdx); 1033 size_t ownerIdx = (itIdx->second)[0]; 1034 int edgeIdxGlo = 0; 1035 size_t faceIdxGlo = mpiRank * nbFaces + nf; 1036 1037 if (myIdx == ownerIdx) 1038 { 1039 CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itIdxGlo = edgeIdx2IdxGlo.find(myIdx); 1040 edgeIdxGlo = (itIdxGlo->second)[0]; 1041 edge_lon(edgeIdxGlo - edge_start) = ( abs( bounds_lon(nv1, nf) - bounds_lon(nv2, nf)) < 180.) ? 1042 (( bounds_lon(nv1, nf) + bounds_lon(nv2, nf)) * 0.5) : 1043 (( bounds_lon(nv1, nf) + bounds_lon(nv2, nf)) * 0.5 -180.); 1044 edge_lat(edgeIdxGlo - edge_start) = ( bounds_lat(nv1, nf) + bounds_lat(nv2, nf) ) * 0.5; 1045 edge_nodes(0, edgeIdxGlo - edge_start) = nodeIdxGlo1; 1046 edge_nodes(1, edgeIdxGlo - edge_start) = nodeIdxGlo2; 1047 } 1074 if (nodeIdxGlo1 != nodeIdxGlo2) 1075 { 1076 CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itIdx = edgeIdx2IdxMin.find(myIdx); 1077 size_t ownerIdx = (itIdx->second)[0]; 1078 int edgeIdxGlo = 0; 1079 size_t faceIdxGlo = nbFacesAccum + nf; 1080 1081 if (myIdx == ownerIdx) 1082 { 1083 CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itIdxGlo = edgeIdx2IdxGlo.find(myIdx); 1084 edgeIdxGlo = (itIdxGlo->second)[0]; 1085 double edgeLon; 1086 double diffLon = abs(bounds_lon(nv1, nf) - bounds_lon(nv2, nf)); 1087 if (diffLon < (180.- prec)) 1088 edgeLon = ( bounds_lon(nv1, nf) + bounds_lon(nv2, nf)) * 0.5; 1089 else if (diffLon > (180.+ prec)) 1090 edgeLon = (bounds_lon(nv1, nf) + bounds_lon(nv2, nf)) * 0.5 -180.; 1091 else 1092 edgeLon = 0.; 1093 edge_lon(edgeIdxGlo - edge_start) = edgeLon; 1094 edge_lat(edgeIdxGlo - edge_start) = ( bounds_lat(nv1, nf) + bounds_lat(nv2, nf) ) * 0.5; 1095 edge_nodes(0, edgeIdxGlo - edge_start) = nodeIdxGlo1; 1096 edge_nodes(1, edgeIdxGlo - edge_start) = nodeIdxGlo2; 1097 } 1098 else 1099 { 1100 CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itIdxGlo = edgeIdxMin2IdxGlo.find(ownerIdx); 1101 edgeIdxGlo = (itIdxGlo->second)[0]; 1102 } 1103 face_edges(nv1,nf) = edgeIdxGlo; 1104 if (edgeIdxGlo2Face.count(edgeIdxGlo) == 0) 1105 { 1106 edgeIdxGloList(iIdx) = edgeIdxGlo; 1107 ++iIdx; 1108 } 1109 edgeIdxGlo2Face[edgeIdxGlo].push_back(faceIdxGlo); 1110 } // nodeIdxGlo1 != nodeIdxGlo2 1048 1111 else 1049 1112 { 1050 CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itIdxGlo = edgeIdxMin2IdxGlo.find(ownerIdx); 1051 edgeIdxGlo = (itIdxGlo->second)[0]; 1052 } 1053 face_edges(nv1,nf) = edgeIdxGlo; 1054 if (edgeIdxGlo2Face.count(edgeIdxGlo) == 0) 1055 { 1056 edgeIdxGloList(iIdx) = edgeIdxGlo; 1057 ++iIdx; 1058 } 1059 edgeIdxGlo2Face[edgeIdxGlo].push_back(faceIdxGlo); 1113 face_edges(nv1,nf) = 999999; 1114 } 1060 1115 } 1061 1116 } … … 1064 1119 // (2.4) Saving remaining variables edge_faces and face_faces 1065 1120 edge_faces.resize(2, edge_count); 1066 face_faces.resize(nvertex, nbFaces );1121 face_faces.resize(nvertex, nbFaces_); 1067 1122 1068 1123 CClientClientDHTSizet dhtEdge2Face (edgeIdxGlo2Face, comm); … … 1071 1126 CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itNodeIdxGlo1, itNodeIdxGlo2; 1072 1127 1073 for (int nf = 0; nf < nbFaces ; ++nf)1128 for (int nf = 0; nf < nbFaces_; ++nf) 1074 1129 { 1075 1130 for (int nv1 = 0; nv1 < nvertex; ++nv1) … … 1097 1152 CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itIdx = edgeIdx2IdxMin.find(myIdx); 1098 1153 size_t ownerIdx = (itIdx->second)[0]; 1099 size_t faceIdxGlo = mpiRank * nbFaces+ nf;1154 size_t faceIdxGlo = nbFacesAccum + nf; 1100 1155 1101 1156 if (myIdx == ownerIdx) … … 1109 1164 edge_faces(0, edgeIdxGlo - edge_start) = face1; 1110 1165 edge_faces(1, edgeIdxGlo - edge_start) = -999; 1111 face_faces(nv1, nf) = -1;1166 face_faces(nv1, nf) = 999999; 1112 1167 } 1113 1168 else … … 1138 1193 vector<size_t> hashValues(4); 1139 1194 CClientClientDHTSizet::Index2VectorInfoTypeMap nodeHash2Idx; 1140 CArray<size_t,1> nodeHashList(nbFaces *nvertex*4);1195 CArray<size_t,1> nodeHashList(nbFaces_*nvertex*4); 1141 1196 size_t iHash = 0; 1142 for (int nf = 0; nf < nbFaces ; ++nf)1197 for (int nf = 0; nf < nbFaces_; ++nf) 1143 1198 { 1144 1199 for (int nv = 0; nv < nvertex; ++nv) … … 1173 1228 CClientClientDHTSizet::Index2VectorInfoTypeMap nodeIdx2IdxMin; 1174 1229 CClientClientDHTSizet::Index2VectorInfoTypeMap nodeIdx2IdxGlo; 1175 CArray<size_t,1> nodeIdxMinList(nbFaces *nvertex*4);1230 CArray<size_t,1> nodeIdxMinList(nbFaces_*nvertex*4); 1176 1231 size_t iIdxMin = 0; 1177 1232 … … 1221 1276 CClientClientDHTSizet::Index2VectorInfoTypeMap edgeHash2Idx; 1222 1277 CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itNodeIdxGlo1, itNodeIdxGlo2; 1223 CArray<size_t,1> edgeHashList(nbFaces*nvertex); 1224 1225 for (int nf = 0; nf < nbFaces; ++nf) 1278 CArray<size_t,1> edgeHashList(nbFaces_*nvertex); 1279 size_t nEdgeHash = 0; 1280 1281 for (int nf = 0; nf < nbFaces_; ++nf) 1226 1282 { 1227 1283 for (int nv1 = 0; nv1 < nvertex; ++nv1) … … 1250 1306 itNodeIdxGlo2 = nodeIdxMin2IdxGlo.find((itNodeIdx2->second)[0]); 1251 1307 nodeIdxGlo2 = (itNodeIdxGlo2->second)[0]; 1252 size_t edgeHash = hashPairOrdered(nodeIdxGlo1, nodeIdxGlo2); 1253 edgeHash2Idx[edgeHash].push_back(generateEdgeIndex(nodeIdxGlo1, nodeIdxGlo2, mpiRank)); 1254 edgeHash2Idx[edgeHash].push_back(mpiRank); 1255 edgeHashList(nf*nvertex + nv1) = edgeHash; 1308 if (nodeIdxGlo1 != nodeIdxGlo2) 1309 { 1310 size_t edgeHash = hashPairOrdered(nodeIdxGlo1, nodeIdxGlo2); 1311 edgeHash2Idx[edgeHash].push_back(generateEdgeIndex(nodeIdxGlo1, nodeIdxGlo2, mpiRank)); 1312 edgeHash2Idx[edgeHash].push_back(mpiRank); 1313 edgeHashList(nEdgeHash) = edgeHash; 1314 ++nEdgeHash; 1315 } 1256 1316 face_nodes(nv1,nf) = nodeIdxGlo1; 1257 1317 } 1258 1318 } 1319 edgeHashList.resizeAndPreserve(nEdgeHash); 1259 1320 1260 1321 // (3.4) Generating global edge indexes … … 1270 1331 CClientClientDHTSizet::Index2VectorInfoTypeMap edgeIdx2IdxMin; 1271 1332 CClientClientDHTSizet::Index2VectorInfoTypeMap edgeIdx2IdxGlo; 1272 CArray<size_t,1> edgeIdxMinList(nbFaces *nvertex);1333 CArray<size_t,1> edgeIdxMinList(nbFaces_*nvertex); 1273 1334 iIdxMin = 0; 1274 1335 … … 1316 1377 edge_lat.resize(edge_count); 1317 1378 edge_nodes.resize(2, edge_count); 1318 face_edges.resize(nvertex, nbFaces );1379 face_edges.resize(nvertex, nbFaces_); 1319 1380 1320 1381 CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator it1, it2; 1321 1382 CClientClientDHTSizet::Index2VectorInfoTypeMap edgeIdxGlo2Face; 1322 CArray<size_t,1> edgeIdxGloList(nbFaces *nvertex);1323 size_t iIdx= 0;1324 1325 for (int nf = 0; nf < nbFaces ; ++nf)1383 CArray<size_t,1> edgeIdxGloList(nbFaces_*nvertex); 1384 size_t nEdge = 0; 1385 1386 for (int nf = 0; nf < nbFaces_; ++nf) 1326 1387 { 1327 1388 for (int nv1 = 0; nv1 < nvertex; ++nv1) … … 1341 1402 size_t nodeIdxGlo2 = (itNodeIdxGlo2->second)[0]; 1342 1403 1343 size_t myIdx = generateEdgeIndex(nodeIdxGlo1, nodeIdxGlo2, mpiRank); 1344 CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itIdx = edgeIdx2IdxMin.find(myIdx); 1345 size_t ownerIdx = (itIdx->second)[0]; 1346 1347 int edgeIdxGlo = 0; 1348 size_t faceIdxGlo = mpiRank * nbFaces + nf; 1349 1350 if (myIdx == ownerIdx) 1351 { 1352 CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itIdxGlo = edgeIdx2IdxGlo.find(myIdx); 1353 edgeIdxGlo = (itIdxGlo->second)[0]; 1354 edge_lon(edgeIdxGlo - edge_start) = ( abs( bounds_lon(nv1, nf) - bounds_lon(nv2, nf)) < 180.) ? 1355 (( bounds_lon(nv1, nf) + bounds_lon(nv2, nf)) * 0.5) : 1356 (( bounds_lon(nv1, nf) + bounds_lon(nv2, nf)) * 0.5 -180.); 1357 edge_lat(edgeIdxGlo-edge_start) = ( bounds_lat(nv1, nf) + bounds_lat(nv2, nf) ) * 0.5; 1358 edge_nodes(0, edgeIdxGlo - edge_start) = nodeIdxGlo1; 1359 edge_nodes(1, edgeIdxGlo - edge_start) = nodeIdxGlo2; 1360 } 1404 if (nodeIdxGlo1 != nodeIdxGlo2) 1405 { 1406 size_t myIdx = generateEdgeIndex(nodeIdxGlo1, nodeIdxGlo2, mpiRank); 1407 CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itIdx = edgeIdx2IdxMin.find(myIdx); 1408 size_t ownerIdx = (itIdx->second)[0]; 1409 int edgeIdxGlo; 1410 size_t faceIdxGlo = nbFacesAccum + nf; 1411 1412 if (myIdx == ownerIdx) 1413 { 1414 CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itIdxGlo = edgeIdx2IdxGlo.find(myIdx); 1415 edgeIdxGlo = (itIdxGlo->second)[0]; 1416 double edgeLon; 1417 double diffLon = abs(bounds_lon(nv1, nf) - bounds_lon(nv2, nf)); 1418 if (diffLon < (180.- prec)) 1419 edgeLon = ( bounds_lon(nv1, nf) + bounds_lon(nv2, nf)) * 0.5; 1420 else if (diffLon > (180.+ prec)) 1421 edgeLon = (bounds_lon(nv1, nf) + bounds_lon(nv2, nf)) * 0.5 -180.; 1422 else 1423 edgeLon = 0.; 1424 edge_lon(edgeIdxGlo - edge_start) = edgeLon; 1425 edge_lat(edgeIdxGlo-edge_start) = ( bounds_lat(nv1, nf) + bounds_lat(nv2, nf) ) * 0.5; 1426 edge_nodes(0, edgeIdxGlo - edge_start) = nodeIdxGlo1; 1427 edge_nodes(1, edgeIdxGlo - edge_start) = nodeIdxGlo2; 1428 } 1429 else 1430 { 1431 CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itIdxGlo = edgeIdxMin2IdxGlo.find(ownerIdx); 1432 edgeIdxGlo = (itIdxGlo->second)[0]; 1433 } 1434 face_edges(nv1,nf) = edgeIdxGlo; 1435 if (edgeIdxGlo2Face.count(edgeIdxGlo) == 0) 1436 { 1437 edgeIdxGloList(nEdge) = edgeIdxGlo; 1438 ++nEdge; 1439 } 1440 edgeIdxGlo2Face[edgeIdxGlo].push_back(faceIdxGlo); 1441 } // nodeIdxGlo1 != nodeIdxGlo2 1361 1442 else 1362 1443 { 1363 CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itIdxGlo = edgeIdxMin2IdxGlo.find(ownerIdx); 1364 edgeIdxGlo = (itIdxGlo->second)[0]; 1365 } 1366 face_edges(nv1,nf) = edgeIdxGlo; 1367 if (edgeIdxGlo2Face.count(edgeIdxGlo) == 0) 1368 { 1369 edgeIdxGloList(iIdx) = edgeIdxGlo; 1370 ++iIdx; 1371 } 1372 edgeIdxGlo2Face[edgeIdxGlo].push_back(faceIdxGlo); 1373 } 1374 } 1375 edgeIdxGloList.resizeAndPreserve(iIdx); 1444 face_edges(nv1,nf) = 999999; 1445 } 1446 } 1447 } 1448 edgeIdxGloList.resizeAndPreserve(nEdge); 1376 1449 1377 1450 // (3.6) Saving remaining variables edge_faces and face_faces 1378 1451 edge_faces.resize(2, edge_count); 1379 face_faces.resize(nvertex, nbFaces );1452 face_faces.resize(nvertex, nbFaces_); 1380 1453 1381 1454 CClientClientDHTSizet dhtEdge2Face (edgeIdxGlo2Face, comm); … … 1383 1456 CClientClientDHTSizet::Index2VectorInfoTypeMap& edgeIdxGlo2FaceIdx = dhtEdge2Face.getInfoIndexMap(); 1384 1457 1385 for (int nf = 0; nf < nbFaces ; ++nf)1458 for (int nf = 0; nf < nbFaces_; ++nf) 1386 1459 { 1387 1460 for (int nv1 = 0; nv1 < nvertex; ++nv1) … … 1394 1467 size_t myNodeIdx1 = generateNodeIndex(hashValues1, mpiRank); 1395 1468 size_t myNodeIdx2 = generateNodeIndex(hashValues2, mpiRank); 1396 it1 = nodeIdx2IdxMin.find(myNodeIdx1); 1397 it2 = nodeIdx2IdxMin.find(myNodeIdx2); 1398 itNodeIdxGlo1 = nodeIdxMin2IdxGlo.find((it1->second)[0]); 1399 itNodeIdxGlo2 = nodeIdxMin2IdxGlo.find((it2->second)[0]); 1400 size_t nodeIdxGlo1 = (itNodeIdxGlo1->second)[0]; 1401 size_t nodeIdxGlo2 = (itNodeIdxGlo2->second)[0]; 1402 1403 size_t myIdx = generateEdgeIndex(nodeIdxGlo1, nodeIdxGlo2, mpiRank); 1404 CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itIdx = edgeIdx2IdxMin.find(myIdx); 1405 size_t ownerIdx = (itIdx->second)[0]; 1406 size_t faceIdxGlo = mpiRank * nbFaces + nf; 1407 1408 if (myIdx == ownerIdx) 1409 { 1410 CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itIdxGlo = edgeIdx2IdxGlo.find(myIdx); 1411 int edgeIdxGlo = (itIdxGlo->second)[0]; 1412 it1 = edgeIdxGlo2FaceIdx.find(edgeIdxGlo); 1413 int face1 = it1->second[0]; 1414 if (it1->second.size() == 1) 1469 if (myNodeIdx1 != myNodeIdx2) 1470 { 1471 it1 = nodeIdx2IdxMin.find(myNodeIdx1); 1472 it2 = nodeIdx2IdxMin.find(myNodeIdx2); 1473 itNodeIdxGlo1 = nodeIdxMin2IdxGlo.find((it1->second)[0]); 1474 itNodeIdxGlo2 = nodeIdxMin2IdxGlo.find((it2->second)[0]); 1475 size_t nodeIdxGlo1 = (itNodeIdxGlo1->second)[0]; 1476 size_t nodeIdxGlo2 = (itNodeIdxGlo2->second)[0]; 1477 1478 size_t myIdx = generateEdgeIndex(nodeIdxGlo1, nodeIdxGlo2, mpiRank); 1479 CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itIdx = edgeIdx2IdxMin.find(myIdx); 1480 size_t ownerIdx = (itIdx->second)[0]; 1481 size_t faceIdxGlo = nbFacesAccum + nf; 1482 1483 if (myIdx == ownerIdx) 1415 1484 { 1416 edge_faces(0, edgeIdxGlo - edge_start) = face1; 1417 edge_faces(1, edgeIdxGlo - edge_start) = -999; 1418 face_faces(nv1, nf) = -1; 1419 } 1485 CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itIdxGlo = edgeIdx2IdxGlo.find(myIdx); 1486 int edgeIdxGlo = (itIdxGlo->second)[0]; 1487 it1 = edgeIdxGlo2FaceIdx.find(edgeIdxGlo); 1488 int face1 = it1->second[0]; 1489 if (it1->second.size() == 1) 1490 { 1491 edge_faces(0, edgeIdxGlo - edge_start) = face1; 1492 edge_faces(1, edgeIdxGlo - edge_start) = -999; 1493 face_faces(nv1, nf) = 999999; 1494 } 1495 else 1496 { 1497 size_t face2 = it1->second[1]; 1498 edge_faces(0, edgeIdxGlo - edge_start) = face1; 1499 edge_faces(1, edgeIdxGlo - edge_start) = face2; 1500 face_faces(nv1, nf) = (faceIdxGlo == face1 ? face2 : face1); 1501 } 1502 } // myIdx == ownerIdx 1420 1503 else 1421 1504 { 1422 size_t face2 = it1->second[1]; 1423 edge_faces(0, edgeIdxGlo - edge_start) = face1; 1424 edge_faces(1, edgeIdxGlo - edge_start) = face2; 1505 CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itIdxGlo = edgeIdxMin2IdxGlo.find(ownerIdx); 1506 size_t edgeIdxGlo = (itIdxGlo->second)[0]; 1507 it1 = edgeIdxGlo2FaceIdx.find(edgeIdxGlo); 1508 int face1 = it1->second[0]; 1509 int face2 = it1->second[1]; 1425 1510 face_faces(nv1, nf) = (faceIdxGlo == face1 ? face2 : face1); 1426 } 1427 } 1511 } // myIdx != ownerIdx 1512 } // myNodeIdx1 != myNodeIdx2 1428 1513 else 1429 { 1430 CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itIdxGlo = edgeIdxMin2IdxGlo.find(ownerIdx); 1431 size_t edgeIdxGlo = (itIdxGlo->second)[0]; 1432 it1 = edgeIdxGlo2FaceIdx.find(edgeIdxGlo); 1433 int face1 = it1->second[0]; 1434 int face2 = it1->second[1]; 1435 face_faces(nv1, nf) = (faceIdxGlo == face1 ? face2 : face1); 1436 } 1514 face_faces(nv1, nf) = 999999; 1437 1515 } 1438 1516 } … … 1443 1521 } // createMeshEpsilon 1444 1522 1523 ///---------------------------------------------------------------- 1524 /*! 1525 * \fn void CMesh::getNghbFacesNodeType(const MPI_Comm& comm, 1526 const CArray<double, 2>& bounds_lon, const CArray<double, 2>& bounds_lat, 1527 CArray<int, 2>& nghbFaces) 1528 * Finds neighboring cells of a local domain for node-type of neighbors. 1529 * \param [in] comm 1530 * \param [in] bounds_lon Array of boundary longitudes. 1531 * \param [in] bounds_lat Array of boundary latitudes. 1532 * \param [out] nghbFaces 2D array of storing global indexes of neighboring cells and their owner procs. 1533 */ 1534 1535 void CMesh::getNghbFacesNodeType(const MPI_Comm& comm, 1536 const CArray<double, 2>& bounds_lon, const CArray<double, 2>& bounds_lat, 1537 CArray<int, 2>& nghbFaces) 1538 { 1539 int nvertex = (bounds_lon.numElements() == 0) ? 1 : bounds_lon.rows(); 1540 int nbFaces = bounds_lon.shape()[1]; 1541 nghbFaces.resize(2, nbFaces*10); // some estimate on max number of neighbouring cells 1542 1543 int mpiRank, mpiSize; 1544 MPI_Comm_rank(comm, &mpiRank); 1545 MPI_Comm_size(comm, &mpiSize); 1546 1547 // For determining the global face index 1548 size_t nbFacesOnProc = nbFaces; 1549 size_t nbFacesAccum; 1550 MPI_Scan(&nbFacesOnProc, &nbFacesAccum, 1, MPI_UNSIGNED_LONG, MPI_SUM, comm); 1551 nbFacesAccum -= nbFaces; 1552 1553 // (1) Generating unique node indexes 1554 // (1.1) Creating a list of hashes for each node and a map nodeHash2Idx <hash, <idx,rank> > 1555 vector<size_t> hashValues(4); 1556 CClientClientDHTSizet::Index2VectorInfoTypeMap nodeHash2Idx; 1557 CArray<size_t,1> nodeHashList(nbFaces*nvertex*4); 1558 size_t iIdx = 0; 1559 for (int nf = 0; nf < nbFaces; ++nf) 1560 { 1561 for (int nv = 0; nv < nvertex; ++nv) 1562 { 1563 hashValues = CMesh::createHashes(bounds_lon(nv, nf), bounds_lat(nv, nf)); 1564 size_t nodeIndex = generateNodeIndex(hashValues, mpiRank); 1565 for (int nh = 0; nh < 4; ++nh) 1566 { 1567 if (nodeHash2Idx.count(hashValues[nh])==0) 1568 { 1569 nodeHash2Idx[hashValues[nh]].push_back(nodeIndex); 1570 nodeHash2Idx[hashValues[nh]].push_back(mpiRank); 1571 nodeHashList(iIdx) = hashValues[nh]; 1572 ++iIdx; 1573 } 1574 } 1575 } 1576 } 1577 nodeHashList.resizeAndPreserve(iIdx); 1578 1579 // (1.2) Generating node indexes 1580 // The ownership criterion: priority of the process holding the smaller index 1581 // Maps generated in this step are: 1582 // nodeHash2Info = <hash, idx1, idx2, idx3....> 1583 // nodeIdx2IdxMin = <idx, idxMin> 1584 // idxMin is a unique node identifier 1585 1586 CClientClientDHTSizet dhtNodeHash(nodeHash2Idx, comm); 1587 dhtNodeHash.computeIndexInfoMapping(nodeHashList); 1588 CClientClientDHTSizet::Index2VectorInfoTypeMap& nodeHash2Info = dhtNodeHash.getInfoIndexMap(); 1589 1590 CClientClientDHTSizet::Index2VectorInfoTypeMap nodeIdx2IdxMin; 1591 1592 for (CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator it = nodeHash2Info.begin(); it != nodeHash2Info.end(); ++it) 1593 { 1594 size_t idxMin = (it->second)[0]; 1595 size_t idx = (it->second)[0]; 1596 for (int i = 2; i < (it->second).size();) 1597 { 1598 if (mpiRank == (it->second)[i+1]) 1599 { 1600 idx = (it->second)[i]; 1601 } 1602 if ((it->second)[i] < idxMin) 1603 { 1604 idxMin = (it->second)[i]; 1605 (it->second)[i] = (it->second)[i-2]; 1606 (it->second)[i+1] = (it->second)[i-1]; 1607 } 1608 i += 2; 1609 } 1610 (it->second)[0] = idxMin; 1611 if (nodeIdx2IdxMin.count(idx) == 0) 1612 { 1613 nodeIdx2IdxMin[idx].push_back(idxMin); 1614 } 1615 } 1616 1617 // (2) Creating maps nodeIdxMin2Face = <nodeIdxMin, [face1, face2, ...]> 1618 CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator it; 1619 CClientClientDHTSizet::Index2VectorInfoTypeMap nodeIdxMin2Face; 1620 CArray<size_t,1> nodeIdxMinList(nbFaces*nvertex*4); 1621 1622 size_t nNode = 0; 1623 1624 for (int nf = 0; nf < nbFaces; ++nf) 1625 { 1626 for (int nv = 0; nv < nvertex; ++nv) 1627 { 1628 vector<size_t> hashValues = CMesh::createHashes(bounds_lon(nv, nf), bounds_lat(nv, nf)); 1629 size_t myNodeIdx = generateNodeIndex(hashValues, mpiRank); 1630 it = nodeIdx2IdxMin.find(myNodeIdx); 1631 size_t nodeIdxMin = (it->second)[0]; 1632 size_t faceIdx = nbFacesAccum + nf; 1633 if (nodeIdxMin2Face.count(nodeIdxMin) == 0) 1634 { 1635 nodeIdxMinList(nNode) = nodeIdxMin; 1636 ++nNode; 1637 } 1638 nodeIdxMin2Face[nodeIdxMin].push_back(faceIdx); 1639 nodeIdxMin2Face[nodeIdxMin].push_back(mpiRank); 1640 } 1641 } 1642 nodeIdxMinList.resizeAndPreserve(nNode); 1643 1644 // (3) Face_face connectivity 1645 1646 // nodeIdxMin2Info = <nodeIdxMin, [face1, face2,...]> 1647 CClientClientDHTSizet dhtNode2Face (nodeIdxMin2Face, comm); 1648 dhtNode2Face.computeIndexInfoMapping(nodeIdxMinList); 1649 CClientClientDHTSizet::Index2VectorInfoTypeMap& nodeIdxMin2Info = dhtNode2Face.getInfoIndexMap(); 1650 CClientClientDHTSizet::Index2VectorInfoTypeMap mapFaces; // auxiliar map 1651 1652 int nbNghb = 0; 1653 CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itNode; 1654 1655 for (int nf = 0; nf < nbFaces; ++nf) 1656 { 1657 for (int nv = 0; nv < nvertex; ++nv) 1658 { 1659 vector<size_t> hashValues = CMesh::createHashes(bounds_lon(nv, nf), bounds_lat(nv, nf)); 1660 size_t myNodeIdx = generateNodeIndex(hashValues, mpiRank); 1661 itNode = nodeIdx2IdxMin.find(myNodeIdx); 1662 size_t nodeIdxMin = (itNode->second)[0]; 1663 1664 itNode = nodeIdxMin2Info.find(nodeIdxMin); 1665 for (int i = 0; i < itNode->second.size();) 1666 { 1667 size_t face = itNode->second[i]; 1668 size_t rank = itNode->second[i+1]; 1669 if (rank != mpiRank) 1670 if (mapFaces.count(face) == 0) 1671 { 1672 nghbFaces(0, nbNghb) = face; 1673 nghbFaces(1, nbNghb) = rank; 1674 ++nbNghb; 1675 mapFaces[face].push_back(face); 1676 } 1677 i += 2; 1678 } 1679 } 1680 } 1681 nghbFaces.resizeAndPreserve(2, nbNghb); 1682 } // getNghbFacesNodeType 1683 1684 ///---------------------------------------------------------------- 1685 /*! 1686 * \fn void CMesh::getNghbFacesEdgeType(const MPI_Comm& comm, 1687 const CArray<double, 2>& bounds_lon, const CArray<double, 2>& bounds_lat, 1688 CArray<int, 2>& nghbFaces) 1689 * Finds neighboring cells of a local domain for edge-type of neighbors. 1690 * \param [in] comm 1691 * \param [in] bounds_lon Array of boundary longitudes. 1692 * \param [in] bounds_lat Array of boundary latitudes. 1693 * \param [out] nghbFaces 2D array of storing global indexes of neighboring cells and their owner procs. 1694 */ 1695 1696 void CMesh::getNghbFacesEdgeType(const MPI_Comm& comm, 1697 const CArray<double, 2>& bounds_lon, const CArray<double, 2>& bounds_lat, 1698 CArray<int, 2>& nghbFaces) 1699 { 1700 int nvertex = (bounds_lon.numElements() == 0) ? 1 : bounds_lon.rows(); 1701 int nbFaces = bounds_lon.shape()[1]; 1702 nghbFaces.resize(2, nbFaces*10); // estimate of max number of neighbouring cells 1703 1704 int mpiRank, mpiSize; 1705 MPI_Comm_rank(comm, &mpiRank); 1706 MPI_Comm_size(comm, &mpiSize); 1707 1708 // For determining the global face index 1709 size_t nbFacesOnProc = nbFaces; 1710 size_t nbFacesAccum; 1711 MPI_Scan(&nbFacesOnProc, &nbFacesAccum, 1, MPI_UNSIGNED_LONG, MPI_SUM, comm); 1712 nbFacesAccum -= nbFaces; 1713 1714 // (1) Generating unique node indexes 1715 // (1.1) Creating a list of hashes for each node and a map nodeHash2Idx <hash, <idx,rank> > 1716 vector<size_t> hashValues(4); 1717 CClientClientDHTSizet::Index2VectorInfoTypeMap nodeHash2Idx; 1718 CArray<size_t,1> nodeHashList(nbFaces*nvertex*4); 1719 size_t iIdx = 0; 1720 for (int nf = 0; nf < nbFaces; ++nf) 1721 { 1722 for (int nv = 0; nv < nvertex; ++nv) 1723 { 1724 hashValues = CMesh::createHashes(bounds_lon(nv, nf), bounds_lat(nv, nf)); 1725 size_t nodeIndex = generateNodeIndex(hashValues, mpiRank); 1726 for (int nh = 0; nh < 4; ++nh) 1727 { 1728 if (nodeHash2Idx.count(hashValues[nh])==0) 1729 { 1730 nodeHash2Idx[hashValues[nh]].push_back(nodeIndex); 1731 nodeHash2Idx[hashValues[nh]].push_back(mpiRank); 1732 nodeHashList(iIdx) = hashValues[nh]; 1733 ++iIdx; 1734 } 1735 } 1736 } 1737 } 1738 nodeHashList.resizeAndPreserve(iIdx); 1739 1740 // (1.2) Generating node indexes 1741 // The ownership criterion: priority of the process holding the smaller index 1742 // Maps generated in this step are: 1743 // nodeHash2Info = <hash, idx1, idx2, idx3....> 1744 // nodeIdx2IdxMin = <idx, idxMin> 1745 // idxMin is a unique node identifier 1746 1747 CClientClientDHTSizet dhtNodeHash(nodeHash2Idx, comm); 1748 dhtNodeHash.computeIndexInfoMapping(nodeHashList); 1749 CClientClientDHTSizet::Index2VectorInfoTypeMap& nodeHash2Info = dhtNodeHash.getInfoIndexMap(); 1750 1751 CClientClientDHTSizet::Index2VectorInfoTypeMap nodeIdx2IdxMin; 1752 1753 for (CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator it = nodeHash2Info.begin(); it != nodeHash2Info.end(); ++it) 1754 { 1755 size_t idxMin = (it->second)[0]; 1756 size_t idx = (it->second)[0]; 1757 for (int i = 2; i < (it->second).size();) 1758 { 1759 if (mpiRank == (it->second)[i+1]) 1760 { 1761 idx = (it->second)[i]; 1762 } 1763 if ((it->second)[i] < idxMin) 1764 { 1765 idxMin = (it->second)[i]; 1766 (it->second)[i] = (it->second)[i-2]; 1767 (it->second)[i+1] = (it->second)[i-1]; 1768 } 1769 i += 2; 1770 } 1771 (it->second)[0] = idxMin; 1772 if (nodeIdx2IdxMin.count(idx) == 0) 1773 { 1774 nodeIdx2IdxMin[idx].push_back(idxMin); 1775 } 1776 } 1777 1778 // (2) Creating map edgeHash2Face = <edgeHash, [[face1, rank1], [face2, rank2]]>, where rank1 = rank2 = ... 1779 1780 CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator it1, it2, it; 1781 CClientClientDHTSizet::Index2VectorInfoTypeMap edgeHash2Face; 1782 CArray<size_t,1> edgeHashList(nbFaces*nvertex); 1783 1784 size_t nEdge = 0; 1785 1786 for (int nf = 0; nf < nbFaces; ++nf) 1787 { 1788 for (int nv1 = 0; nv1 < nvertex; ++nv1) 1789 { 1790 // Getting indexes of edge's nodes 1791 int nv2 = (nv1 < nvertex -1 ) ? (nv1 + 1) : (nv1 + 1 - nvertex); // cyclic rotation 1792 vector<size_t> hashValues1 = CMesh::createHashes(bounds_lon(nv1, nf), bounds_lat(nv1, nf)); 1793 vector<size_t> hashValues2 = CMesh::createHashes(bounds_lon(nv2, nf), bounds_lat(nv2, nf)); 1794 size_t myNodeIdx1 = generateNodeIndex(hashValues1, mpiRank); 1795 size_t myNodeIdx2 = generateNodeIndex(hashValues2, mpiRank); 1796 it1 = nodeIdx2IdxMin.find(myNodeIdx1); 1797 it2 = nodeIdx2IdxMin.find(myNodeIdx2); 1798 size_t nodeIdxMin1 = (it1->second)[0]; 1799 size_t nodeIdxMin2 = (it2->second)[0]; 1800 size_t faceIdx = nbFacesAccum + nf; 1801 1802 if (nodeIdxMin1 != nodeIdxMin2) 1803 { 1804 size_t edgeHash = hashPairOrdered(nodeIdxMin1, nodeIdxMin2); 1805 if (edgeHash2Face.count(edgeHash) == 0) 1806 { 1807 edgeHashList(nEdge) = edgeHash; 1808 ++nEdge; 1809 } 1810 edgeHash2Face[edgeHash].push_back(faceIdx); 1811 edgeHash2Face[edgeHash].push_back(mpiRank); 1812 } // nodeIdxMin1 != nodeIdxMin2 1813 } 1814 } 1815 edgeHashList.resizeAndPreserve(nEdge); 1816 1817 // (3) Face_face connectivity 1818 1819 int nbNghb = 0; 1820 CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itNode1, itNode2; 1821 1822 // edgeHash2Info = <edgeHash, [[face1, rank1], [face2, rank2]]> 1823 CClientClientDHTSizet dhtEdge2Face (edgeHash2Face, comm); 1824 dhtEdge2Face.computeIndexInfoMapping(edgeHashList); 1825 CClientClientDHTSizet::Index2VectorInfoTypeMap& edgeHash2Info = dhtEdge2Face.getInfoIndexMap(); 1826 CClientClientDHTSizet::Index2VectorInfoTypeMap mapFaces; // auxiliar map 1827 1828 for (int nf = 0; nf < nbFaces; ++nf) 1829 { 1830 nbNghb = 0; 1831 for (int nv1 = 0; nv1 < nvertex; ++nv1) 1832 { 1833 // Getting indexes of edge's nodes 1834 int nv2 = (nv1 < nvertex -1 ) ? (nv1 + 1) : (nv1 + 1 - nvertex); // cyclic rotation 1835 vector<size_t> hashValues1 = CMesh::createHashes(bounds_lon(nv1, nf), bounds_lat(nv1, nf)); 1836 vector<size_t> hashValues2 = CMesh::createHashes(bounds_lon(nv2, nf), bounds_lat(nv2, nf)); 1837 1838 size_t myNodeIdx1 = generateNodeIndex(hashValues1, mpiRank); 1839 size_t myNodeIdx2 = generateNodeIndex(hashValues2, mpiRank); 1840 itNode1 = nodeIdx2IdxMin.find(myNodeIdx1); 1841 itNode2 = nodeIdx2IdxMin.find(myNodeIdx2); 1842 size_t nodeIdxMin1 = (itNode1->second)[0]; 1843 size_t nodeIdxMin2 = (itNode2->second)[0]; 1844 1845 if (nodeIdxMin1 != nodeIdxMin2) 1846 { 1847 size_t edgeHash = hashPairOrdered(nodeIdxMin1, nodeIdxMin2); 1848 it1 = edgeHash2Info.find(edgeHash); 1849 1850 for (int i = 0; i < it1->second.size();) 1851 { 1852 size_t face = it1->second[i]; 1853 size_t rank = it1->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 } // nodeIdxMin1 != nodeIdxMin2 1865 } 1866 } 1867 nghbFaces.resizeAndPreserve(2, nbNghb); 1868 } // getNghbFacesEdgeType 1869 1870 ///---------------------------------------------------------------- 1871 /*! 1872 * \fn void getDomainExtended(const int nghbType, const MPI_Comm& comm, 1873 const CArray<double, 2>& bounds_lon, const CArray<double, 2>& bounds_lat, 1874 CArray<size_t, 1>& nghbFaces) 1875 * \param [in] nghbType 0 for faces sharing nodes, otherwise for faces sharing edges. 1876 * \param [in] comm 1877 * \param [in] bounds_lon Array of boundary longitudes. 1878 * \param [in] bounds_lat Array of boundary latitudes. 1879 * \param [out] nghbFaces Array containing neighboring faces and their ranks. The ordering is face1, rank1,.... 1880 */ 1881 1882 void CMesh::getDomainExtended(const int nghbType, const MPI_Comm& comm, 1883 const CArray<double, 2>& bounds_lon, const CArray<double, 2>& bounds_lat, 1884 CArray<int, 2>& nghbFaces) 1885 { 1886 if (nghbType == 0) 1887 getNghbFacesNodeType(comm, bounds_lon, bounds_lat, nghbFaces); 1888 else 1889 getNghbFacesEdgeType(comm, bounds_lon, bounds_lat, nghbFaces); 1890 } // getDomainExtended 1891 1445 1892 } // namespace xios -
XIOS/trunk/src/node/mesh.hpp
r924 r929 57 57 CArray<int, 2> face_faces; 58 58 59 void createMesh(const CArray<double, 1>&, const CArray<double, 1>&, 59 void createMesh(const CArray<double, 1>&, const CArray<double, 1>&, 60 60 const CArray<double, 2>&, const CArray<double, 2>& ); 61 61 … … 63 63 const CArray<double, 1>&, const CArray<double, 1>&, 64 64 const CArray<double, 2>&, const CArray<double, 2>& ); 65 66 void getDomainExtended(const int, const MPI_Comm&, 67 const CArray<double, 2>&, const CArray<double, 2>&, 68 CArray<int, 2>&); 65 69 66 70 static CMesh* getMesh(StdString, int); … … 68 72 private: 69 73 70 int nbNodes ;71 int nbEdges ;72 int nbFaces ;74 int nbNodes_; 75 int nbEdges_; 76 int nbFaces_; 73 77 74 78 static std::map <StdString, CMesh> meshList; … … 76 80 CClientClientDHTSizet* pNodeGlobalIndex; // pointer to a map <nodeHash, nodeIdxGlo> 77 81 CClientClientDHTSizet* pEdgeGlobalIndex; // pointer to a map <edgeHash, edgeIdxGlo> 82 void getNghbFacesNodeType(const MPI_Comm&, const CArray<double, 2>&, const CArray<double, 2>&, CArray<int, 2>&); 83 void getNghbFacesEdgeType(const MPI_Comm&, const CArray<double, 2>&, const CArray<double, 2>&, CArray<int, 2>&); 78 84 79 85 vector<size_t> createHashes (const double, const double);
Note: See TracChangeset
for help on using the changeset viewer.