546 ierr = MPI_Comm_size(PROTEUS_COMM_WORLD,&size);
547 ierr = MPI_Comm_rank(PROTEUS_COMM_WORLD,&rank);
577 valarray<int> nodeOffsets_old(size+1);
578 nodeOffsets_old[0] = 0;
579 for (
int sdN=0; sdN < size; sdN++)
581 nodeOffsets_old[sdN+1] = nodeOffsets_old[sdN] +
588 int nNodes_subdomain = (nodeOffsets_old[rank+1] - nodeOffsets_old[rank]);
590 PetscInt *nodeNeighborsOffsets_subdomain,*nodeNeighbors_subdomain,*weights_subdomain;
591 PetscMalloc(
sizeof(PetscInt)*(nNodes_subdomain+1),&nodeNeighborsOffsets_subdomain);
594 nodeNeighborsOffsets_subdomain[0] = 0;
595 for (
int nN = 0,offset=0; nN < nNodes_subdomain; nN++)
597 int nN_global = nodeOffsets_old[rank] + nN;
602 nodeNeighbors_subdomain[offset++] = mesh.
nodeStarArray[offset_global];
605 nodeNeighborsOffsets_subdomain[nN+1]=offset;
606 sort(&nodeNeighbors_subdomain[nodeNeighborsOffsets_subdomain[nN]],&nodeNeighbors_subdomain[nodeNeighborsOffsets_subdomain[nN+1]]);
607 int weight= (nodeNeighborsOffsets_subdomain[nN+1] - nodeNeighborsOffsets_subdomain[nN]);
608 for (
int k=nodeNeighborsOffsets_subdomain[nN];k<nodeNeighborsOffsets_subdomain[nN+1];k++)
609 weights_subdomain[k] = weight;
621 ierr = MatCreateMPIAdj(PROTEUS_COMM_WORLD,
624 nodeNeighborsOffsets_subdomain,
625 nodeNeighbors_subdomain,
627 &petscAdjacency);CHKERRABORT(PROTEUS_COMM_WORLD, ierr);
628 MatPartitioning petscPartition;
629 MatPartitioningCreate(PROTEUS_COMM_WORLD,&petscPartition);
630 MatPartitioningSetAdjacency(petscPartition,petscAdjacency);
631 MatPartitioningSetFromOptions(petscPartition);
634 IS nodePartitioningIS_new;
635 MatPartitioningApply(petscPartition,&nodePartitioningIS_new);
636 MatPartitioningDestroy(&petscPartition);
639 valarray<int> nNodes_subdomain_new(size);
640 ISPartitioningCount(nodePartitioningIS_new,size,&nNodes_subdomain_new[0]);
643 valarray<int> nodeOffsets_new(size+1);
644 nodeOffsets_new[0] = 0;
645 for (
int sdN = 0; sdN < size; sdN++)
646 nodeOffsets_new[sdN+1] = nodeOffsets_new[sdN] + nNodes_subdomain_new[sdN];
649 IS nodeNumberingIS_subdomain_old2new;
650 ISPartitioningToNumbering(nodePartitioningIS_new,&nodeNumberingIS_subdomain_old2new);
655 IS nodeNumberingIS_global_old2new;
656 ISAllGather(nodeNumberingIS_subdomain_old2new,&nodeNumberingIS_global_old2new);
658 const PetscInt * nodeNumbering_global_old2new;
659 ISGetIndices(nodeNumberingIS_global_old2new,&nodeNumbering_global_old2new);
662 valarray<int> nodeNumbering_global_new2old(mesh.
nNodes_global);
664 nodeNumbering_global_new2old[nodeNumbering_global_old2new[nN]] = nN;
699 MPI_Recv(elementMask,PetscBTLength(mesh.
nElements_global),MPI_CHAR,rank-1,0,PROTEUS_COMM_WORLD,&status);
702 set<int> elements_subdomain_owned;
703 for (
int nN = nodeOffsets_new[rank]; nN < nodeOffsets_new[rank+1]; nN++)
705 int nN_global_old = nodeNumbering_global_new2old[nN];
710 if (!PetscBTLookupSet(elementMask,eN_star_old))
712 elements_subdomain_owned.insert(eN_star_old);
743 MPI_Send(elementMask,PetscBTLength(mesh.
nElements_global),MPI_CHAR,rank+1,0,PROTEUS_COMM_WORLD);
744 ierr = PetscBTDestroy(&elementMask);
746 cerr<<
"Error in PetscBTDestroy"<<endl;
751 valarray<int> nElements_subdomain_new(size),
752 elementOffsets_new(size+1);
753 for (
int sdN = 0; sdN < size; sdN++)
756 nElements_subdomain_new[sdN] = int(elements_subdomain_owned.size());
758 nElements_subdomain_new[sdN] = 0;
760 valarray<int> nElements_subdomain_new_send = nElements_subdomain_new;
761 MPI_Allreduce(&nElements_subdomain_new_send[0],&nElements_subdomain_new[0],size,MPI_INT,MPI_SUM,PROTEUS_COMM_WORLD);
763 elementOffsets_new[0] = 0;
764 for (
int sdN = 0; sdN < size; sdN++)
765 elementOffsets_new[sdN+1] = elementOffsets_new[sdN] + nElements_subdomain_new[sdN];
768 valarray<int> elementNumbering_subdomain_new2old(elements_subdomain_owned.size());
769 set<int>::iterator eN_ownedp = elements_subdomain_owned.begin();
770 for (
int eN = 0; eN < int(elements_subdomain_owned.size()); eN++)
772 elementNumbering_subdomain_new2old[eN] = *eN_ownedp++;
775 IS elementNumberingIS_subdomain_new2old;
776 ISCreateGeneral(PROTEUS_COMM_WORLD,elements_subdomain_owned.size(),&elementNumbering_subdomain_new2old[0],PETSC_COPY_VALUES,
777 &elementNumberingIS_subdomain_new2old);
778 IS elementNumberingIS_global_new2old;
779 ISAllGather(elementNumberingIS_subdomain_new2old,&elementNumberingIS_global_new2old);
781 const PetscInt *elementNumbering_global_new2old;
782 ISGetIndices(elementNumberingIS_global_new2old,&elementNumbering_global_new2old);
787 elementNumbering_global_old2new[elementNumbering_global_new2old[eN]] = eN;
802 MPI_Status status_elementBoundaries;
803 PetscBT elementBoundaryMask;
807 MPI_Recv(elementBoundaryMask,PetscBTLength(mesh.
nElementBoundaries_global),MPI_CHAR,rank-1,0,PROTEUS_COMM_WORLD,&status_elementBoundaries);
811 set<int> elementBoundaries_subdomain_owned;
815 int lface[6][4] = {{0,1,2,3},
822 for (
int nN = nodeOffsets_new[rank]; nN < nodeOffsets_new[rank+1]; nN++)
824 int nN_global_old = nodeNumbering_global_new2old[nN];
830 int eN_star_new = elementNumbering_global_old2new[eN_star_old];
836 bool foundNode =
false;
840 if (nN_global_old_across == nN_global_old) foundNode =
true;
845 if (!PetscBTLookupSet(elementBoundaryMask,ebN_global))
846 elementBoundaries_subdomain_owned.insert(ebN_global);
856 for (
int nN = nodeOffsets_new[rank]; nN < nodeOffsets_new[rank+1]; nN++)
858 int nN_global_old = nodeNumbering_global_new2old[nN];
864 int eN_star_new = elementNumbering_global_old2new[eN_star_old];
871 if (nN_global_old_across != nN_global_old)
874 if (!PetscBTLookupSet(elementBoundaryMask,ebN_global))
875 elementBoundaries_subdomain_owned.insert(ebN_global);
885 ierr = PetscBTDestroy(&elementBoundaryMask);
887 cerr<<
"Error in PetscBTDestroy for elementBoundaries"<<endl;
889 valarray<int> nElementBoundaries_subdomain_new(size),
890 elementBoundaryOffsets_new(size+1);
891 for (
int sdN=0;sdN<size;sdN++)
893 nElementBoundaries_subdomain_new[sdN] = elementBoundaries_subdomain_owned.size();
895 nElementBoundaries_subdomain_new[sdN] = 0;
896 valarray<int> nElementBoundaries_subdomain_new_send=nElementBoundaries_subdomain_new;
897 MPI_Allreduce(&nElementBoundaries_subdomain_new_send[0],&nElementBoundaries_subdomain_new[0],size,MPI_INT,MPI_SUM,PROTEUS_COMM_WORLD);
898 elementBoundaryOffsets_new[0] = 0;
899 for (
int sdN=0;sdN<size;sdN++)
900 elementBoundaryOffsets_new[sdN+1] = elementBoundaryOffsets_new[sdN]+nElementBoundaries_subdomain_new[sdN];
905 valarray<int> elementBoundaryNumbering_new2old(elementBoundaries_subdomain_owned.size());
906 set<int>::iterator ebN_ownedp=elementBoundaries_subdomain_owned.begin();
907 for (
int ebN=0;ebN<int(elementBoundaries_subdomain_owned.size());ebN++)
909 elementBoundaryNumbering_new2old[ebN] = *ebN_ownedp++;
911 IS elementBoundaryNumberingIS_subdomain_new2old;
912 ISCreateGeneral(PROTEUS_COMM_WORLD,elementBoundaries_subdomain_owned.size(),&elementBoundaryNumbering_new2old[0],PETSC_COPY_VALUES,&elementBoundaryNumberingIS_subdomain_new2old);
913 IS elementBoundaryNumberingIS_global_new2old;
914 ISAllGather(elementBoundaryNumberingIS_subdomain_new2old,&elementBoundaryNumberingIS_global_new2old);
915 const PetscInt *elementBoundaryNumbering_global_new2old;
917 ISGetIndices(elementBoundaryNumberingIS_global_new2old,&elementBoundaryNumbering_global_new2old);
920 elementBoundaryNumbering_old2new_global[elementBoundaryNumbering_global_new2old[ebN]] = ebN;
937 map<NodeTuple<2>,
int> nodesEdgeMap_global;
938 set<int> edges_subdomain_owned;
943 const int nN0_global = nodeNumbering_global_old2new[nN0_global_old];
944 const int nN1_global = nodeNumbering_global_old2new[nN1_global_old];
946 nodes[0] = nN0_global;
947 nodes[1] = nN1_global;
949 nodesEdgeMap_global[et] = ig;
950 if (nodeOffsets_new[rank] <= et.
nodes[0] && et.
nodes[0] < nodeOffsets_new[rank+1])
951 edges_subdomain_owned.insert(ig);
953 valarray<int> nEdges_subdomain_new(size),
954 edgeOffsets_new(size+1);
956 for (
int sdN=0; sdN < size; sdN++)
958 nEdges_subdomain_new[sdN] = edges_subdomain_owned.size();
960 nEdges_subdomain_new[sdN] = 0;
962 valarray<int> nEdges_subdomain_new_send=nEdges_subdomain_new;
963 MPI_Allreduce(&nEdges_subdomain_new_send[0],&nEdges_subdomain_new[0],size,MPI_INT,MPI_SUM,PROTEUS_COMM_WORLD);
964 edgeOffsets_new[0] = 0;
965 for (
int sdN=0;sdN<size;sdN++)
966 edgeOffsets_new[sdN+1] = edgeOffsets_new[sdN]+nEdges_subdomain_new[sdN];
969 valarray<int> edgeNumbering_new2old(edges_subdomain_owned.size());
970 set<int>::iterator edges_ownedp = edges_subdomain_owned.begin();
971 for (
int i=0; i < int(edges_subdomain_owned.size());i++)
972 edgeNumbering_new2old[i] = *edges_ownedp++;
974 IS edgeNumberingIS_subdomain_new2old;
975 ISCreateGeneral(PROTEUS_COMM_WORLD,edges_subdomain_owned.size(),&edgeNumbering_new2old[0],PETSC_COPY_VALUES,&edgeNumberingIS_subdomain_new2old);
976 IS edgeNumberingIS_global_new2old;
977 ISAllGather(edgeNumberingIS_subdomain_new2old,&edgeNumberingIS_global_new2old);
978 const PetscInt *edgeNumbering_global_new2old;
979 valarray<int> edgeNumbering_old2new_global(mesh.
nEdges_global);
980 ISGetIndices(edgeNumberingIS_global_new2old,&edgeNumbering_global_new2old);
983 edgeNumbering_old2new_global[edgeNumbering_global_new2old[ig]] = ig;
988 valarray<int> edgeNodesArray_newNodesAndEdges(2*mesh.
nEdges_global);
989 map<NodeTuple<2>,
int > nodesEdgeMap_global_new;
994 const int nN0_global = nodeNumbering_global_old2new[nN0_global_old];
995 const int nN1_global = nodeNumbering_global_old2new[nN1_global_old];
997 const int edge_new = edgeNumbering_old2new_global[ig];
998 edgeNodesArray_newNodesAndEdges[edge_new*2+0] = nN0_global;
999 edgeNodesArray_newNodesAndEdges[edge_new*2+1] = nN1_global;
1001 nodes[0] = nN0_global;
1002 nodes[1] = nN1_global;
1004 nodesEdgeMap_global_new[et] = edge_new;
1015 set<int> elements_overlap,nodes_overlap,elementBoundaries_overlap,edges_overlap;
1016 for (
int nN = nodeOffsets_new[rank]; nN < nodeOffsets_new[rank+1]; nN++)
1018 int nN_global_old = nodeNumbering_global_new2old[nN];
1024 int nN_neig_new = nodeNumbering_global_old2new[nN_neig_old];
1026 bool offproc = nN_neig_new >= nodeOffsets_new[rank+1] || nN_neig_new < nodeOffsets_new[rank];
1028 nodes_overlap.insert(nN_neig_new);
1035 int eN_star_new = elementNumbering_global_old2new[eN_star_old];
1036 bool offproc = eN_star_new >= elementOffsets_new[rank+1] || eN_star_new < elementOffsets_new[rank];
1038 elements_overlap.insert(eN_star_new);
1046 if (ebN_global < elementBoundaryOffsets_new[rank] || ebN_global >= elementBoundaryOffsets_new[rank+1])
1048 elementBoundaries_overlap.insert(ebN_global);
1056 const int nN0_global = nodeNumbering_global_old2new[nN0_global_old];
1057 const int nN1_global = nodeNumbering_global_old2new[nN1_global_old];
1058 bool foundEdge =
false;
1060 nodes[0] = nN0_global;
1061 nodes[1] = nN1_global;
1063 const int edge_global = nodesEdgeMap_global_new[et];
1064 if (edge_global < edgeOffsets_new[rank] || edge_global >= edgeOffsets_new[rank+1])
1065 edges_overlap.insert(edge_global);
1074 int overlap_remaining = nNodes_overlap -1;
1076 set<int> last_nodes_added2overlap = nodes_overlap;
1077 while (overlap_remaining > 0)
1079 set<int> new_nodes_overlap,new_elements_overlap,new_elementBoundaries_overlap,new_edges_overlap;
1080 set<int>::iterator nN_p = last_nodes_added2overlap.begin();
1081 while (nN_p != last_nodes_added2overlap.end())
1083 int nN_global_new = *nN_p;
1084 int nN_global_old = nodeNumbering_global_new2old[nN_global_new];
1089 int nN_neig_new = nodeNumbering_global_old2new[nN_neig_old];
1092 bool offproc = nN_neig_new >= nodeOffsets_new[rank+1] || nN_neig_new < nodeOffsets_new[rank];
1094 new_nodes_overlap.insert(nN_neig_new);
1100 set<int>::iterator nN_newp = last_nodes_added2overlap.begin();
1101 while (nN_newp != last_nodes_added2overlap.end())
1103 int nN_global_new = *nN_newp;
1104 int nN_global_old = nodeNumbering_global_new2old[nN_global_new];
1109 int eN_star_new = elementNumbering_global_old2new[eN_star_old];
1110 bool offproc = eN_star_new >= elementOffsets_new[rank+1] || eN_star_new < elementOffsets_new[rank];
1112 new_elements_overlap.insert(eN_star_new);
1121 if (ebN_global < elementBoundaryOffsets_new[rank] || ebN_global >= elementBoundaryOffsets_new[rank+1])
1122 new_elementBoundaries_overlap.insert(ebN_global);
1130 const int nN0_global = nodeNumbering_global_old2new[nN0_global_old];
1131 const int nN1_global = nodeNumbering_global_old2new[nN1_global_old];
1132 bool foundEdge =
false;
1134 nodes[0] = nN0_global;
1135 nodes[1] = nN1_global;
1137 const int edge_global = nodesEdgeMap_global_new[et];
1138 if (edge_global < edgeOffsets_new[rank] || edge_global >= edgeOffsets_new[rank+1])
1139 new_edges_overlap.insert(edge_global);
1144 last_nodes_added2overlap.clear();
1145 set_difference(new_nodes_overlap.begin(),new_nodes_overlap.end(),
1146 nodes_overlap.begin(),nodes_overlap.end(),
1147 std::inserter(last_nodes_added2overlap,last_nodes_added2overlap.begin()));
1150 for (set<int>::iterator nN_addedp = new_nodes_overlap.begin();
1151 nN_addedp != new_nodes_overlap.end();
1154 nodes_overlap.insert(*nN_addedp);
1156 for (set<int>::iterator eN_addedp = new_elements_overlap.begin();
1157 eN_addedp != new_elements_overlap.end();
1160 elements_overlap.insert(*eN_addedp);
1162 for (set<int>::iterator ebN_addedp = new_elementBoundaries_overlap.begin();
1163 ebN_addedp != new_elementBoundaries_overlap.end();
1166 elementBoundaries_overlap.insert(*ebN_addedp);
1168 for (set<int>::iterator edge_addedp = new_edges_overlap.begin();
1169 edge_addedp != new_edges_overlap.end();
1172 edges_overlap.insert(*edge_addedp);
1176 overlap_remaining--;
1196 map<int,int> nodeNumbering_global2subdomain;
1197 map<int,int> elementBoundaryNumbering_global2subdomain;
1202 for (
int nN = 0; nN < nNodes_subdomain_new[rank]; nN++)
1204 int nN_global_new = nN + nodeOffsets_new[rank];
1205 int nN_global_old = nodeNumbering_global_new2old[nN_global_new];
1206 nodeNumbering_subdomain2global[nN] = nN_global_new;
1207 nodeNumbering_global2subdomain[nN_global_new] = nN;
1224 set<int>::iterator nN_p = nodes_overlap.begin();
1225 for (
int nN = nNodes_subdomain_new[rank]; nN < nNodes_subdomain_new[rank] + int(nodes_overlap.size()); nN++)
1227 int nN_global_new = *nN_p++;
1228 int nN_global_old = nodeNumbering_global_new2old[nN_global_new];
1229 nodeNumbering_subdomain2global[nN] = nN_global_new;
1230 nodeNumbering_global2subdomain[nN_global_new] = nN;
1239 for (
int eN = 0; eN < nElements_subdomain_new[rank]; eN++)
1241 int eN_global_new = elementOffsets_new[rank] + eN;
1242 int eN_global_old = elementNumbering_global_new2old[eN_global_new];
1243 elementNumbering_subdomain2global[eN] = eN_global_new;
1248 int nN_global_new = nodeNumbering_global_old2new[nN_global_old];
1249 int nN_subdomain = nodeNumbering_global2subdomain[nN_global_new];
1255 set<int>::iterator eN_p = elements_overlap.begin();
1256 for (
int eN = nElements_subdomain_new[rank]; eN < nElements_subdomain_new[rank] + int(elements_overlap.size()); eN++)
1258 int eN_global_new = *eN_p++;
1259 int eN_global_old = elementNumbering_global_new2old[eN_global_new];
1260 elementNumbering_subdomain2global[eN] = eN_global_new;
1265 int nN_global_new = nodeNumbering_global_old2new[nN_global_old];
1266 int nN_subdomain = nodeNumbering_global2subdomain[nN_global_new];
1273 for (
int ebN=0; ebN < nElementBoundaries_subdomain_new[rank]; ebN++)
1275 int ebN_global = ebN + elementBoundaryOffsets_new[rank];
1276 elementBoundaryNumbering_subdomain2global[ebN]=ebN_global;
1277 elementBoundaryNumbering_global2subdomain[ebN_global] = ebN;
1280 set<int>::iterator ebN_p = elementBoundaries_overlap.begin();
1281 for(
int ebN=nElementBoundaries_subdomain_new[rank];ebN < nElementBoundaries_subdomain_new[rank] + int(elementBoundaries_overlap.size()); ebN++)
1283 int ebN_global = *ebN_p++;
1284 elementBoundaryNumbering_subdomain2global[ebN] = ebN_global;
1285 elementBoundaryNumbering_global2subdomain[ebN_global] = ebN;
1290 for (
int eN=0;eN<nElements_subdomain_new[rank];eN++)
1292 int eN_global = eN+elementOffsets_new[rank];
1298set<int>::iterator eN_p2 = elements_overlap.begin();
1299for (
int eN = nElements_subdomain_new[rank]; eN < nElements_subdomain_new[rank] + int(elements_overlap.size()); eN++)
1301 int eN_global_new = *eN_p2++;
1315for (
int i=0; i < nEdges_subdomain_new[rank]; i++)
1317 const int ig = i+edgeOffsets_new[rank];
1318 const int nN0_global = edgeNodesArray_newNodesAndEdges[ig*2+0];
1319 const int nN1_global = edgeNodesArray_newNodesAndEdges[ig*2+1];
1321 const int nN0_subdomain = nodeNumbering_global2subdomain[nN0_global];
1322 const int nN1_subdomain = nodeNumbering_global2subdomain[nN1_global];
1325 edgeNumbering_subdomain2global[i] = ig;
1328set<int>::iterator edge_p = edges_overlap.begin();
1329for (
int i=nEdges_subdomain_new[rank]; i < nEdges_subdomain_new[rank] + int(edges_overlap.size()); i++)
1331 const int ig =*edge_p++;
1332 const int nN0_global = edgeNodesArray_newNodesAndEdges[ig*2+0];
1333 const int nN1_global = edgeNodesArray_newNodesAndEdges[ig*2+1];
1335 const int nN0_subdomain = nodeNumbering_global2subdomain[nN0_global];
1336 const int nN1_subdomain = nodeNumbering_global2subdomain[nN1_global];
1339 edgeNumbering_subdomain2global[i] = ig;
1405 int eN_global_new = elementNumbering_subdomain2global[eN];
1406 int eN_global_old = elementNumbering_global_new2old[eN_global_new];
1433 for (
int sdN = 0; sdN < size+1; sdN++)
1464 ISRestoreIndices(nodeNumberingIS_global_old2new,&nodeNumbering_global_old2new);
1466 ISDestroy(&nodePartitioningIS_new);
1467 ISDestroy(&nodeNumberingIS_subdomain_old2new);
1468 ISDestroy(&nodeNumberingIS_global_old2new);
1470 ISRestoreIndices(elementNumberingIS_global_new2old,&elementNumbering_global_new2old);
1472 ISDestroy(&elementNumberingIS_subdomain_new2old);
1473 ISDestroy(&elementNumberingIS_global_new2old);
1475 ISRestoreIndices(elementBoundaryNumberingIS_global_new2old,&elementBoundaryNumbering_global_new2old);
1477 ISDestroy(&elementBoundaryNumberingIS_subdomain_new2old);
1478 ISDestroy(&elementBoundaryNumberingIS_global_new2old);
1480 ISRestoreIndices(edgeNumberingIS_global_new2old,&edgeNumbering_global_new2old);
1482 ISDestroy(&edgeNumberingIS_subdomain_new2old);
1483 ISDestroy(&edgeNumberingIS_global_new2old);
1489 Mesh& newMesh,
int nNodes_overlap,
double memHardLimit)
1491 using namespace std;
1492 PetscErrorCode ierr;
1493 PetscMPIInt size,rank;
1495 ierr = MPI_Comm_size(PROTEUS_COMM_WORLD,&size);CHKERRABORT(PROTEUS_COMM_WORLD, ierr);
1496 ierr = MPI_Comm_rank(PROTEUS_COMM_WORLD,&rank);CHKERRABORT(PROTEUS_COMM_WORLD, ierr);
1498 PetscLogStage partitioning_stage;
1499 PetscLogStageRegister(
"Mesh Partition",&partitioning_stage);
1500 PetscLogStagePush(partitioning_stage);
1527 bool failed =
false;
1528 const int simplexDim = 4;
1529 const int vertexDim = 3;
1531 std::string vertexFileName = std::string(filebase) +
".node" ;
1532 std::string elementFileName = std::string(filebase) +
".ele" ;
1533 std::string elementBoundaryFileName = std::string(filebase) +
".face" ;
1534 std::string edgeFileName = std::string(filebase) +
".edge" ;
1544 logEvent(
"Building default partitioning",5);
1545 int read_elements_event;
1546 PetscLogEventRegister(
"Read eles",0,&read_elements_event);
1547 PetscLogEventBegin(read_elements_event,0,0,0,0);
1548 std::ifstream vertexFile(vertexFileName.c_str());
1549 if (!vertexFile.good())
1551 std::cerr<<
"cannot open Tetgen node file "
1552 <<vertexFileName<<std::endl;
1556 int hasVertexMarkers(0),hasVertexAttributes(0),nSpace(3),nNodes_global;
1558 vertexFile >>
eatcomments >> nNodes_global >> nSpace >> hasVertexAttributes >> hasVertexMarkers >>
eatline ;
1559 assert(nNodes_global > 0);
1560 assert(nSpace == 3);
1566 if (hasVertexAttributes > 0 && rank==0)
1568 std::cerr<<
"WARNING Tetgen nodes hasAttributes= "<<hasVertexAttributes
1569 <<
" > 0 will treat first value as integer id for boundary!!"<<std::endl;
1570 hasVertexMarkers = 1;
1578 valarray<int> nodeOffsets_old(size+1);
1579 nodeOffsets_old[0] = 0;
1580 for (
int sdN=0; sdN < size; sdN++)
1582 nodeOffsets_old[sdN+1] = nodeOffsets_old[sdN] +
1583 int(nNodes_global)/size + (int(nNodes_global)%size > sdN);
1585 int nNodes_subdomain_old = nodeOffsets_old[rank+1] - nodeOffsets_old[rank];
1586 int nNodes_subdomain_max = int(nNodes_global)/size + (int(nNodes_global)%size > 0);
1587 char log_buffer[100];
1588 sprintf(log_buffer,
"First partitioning: Max nNodes_subdomain %d nNodes_global %d",nNodes_subdomain_max,nNodes_global);
1589 logEvent(log_buffer,5);
1595 logEvent(
"Determining nodal connectivity to construct new partitioning",5);
1596 std::ifstream elementFile(elementFileName.c_str());
1597 if (!elementFile.good())
1599 std::cerr<<
"cannot open Tetgen element file "
1600 <<elementFileName<<std::endl;
1605 int nNodesPerSimplex(simplexDim),hasElementMarkers = 0,nElements_global;
1606 elementFile >>
eatcomments >> nElements_global >> nNodesPerSimplex >> hasElementMarkers >>
eatline;
1607 assert(nElements_global > 0);
1608 assert(nNodesPerSimplex == simplexDim);
1610 valarray<int> element_nodes_old(4);
1612 map<int,valarray<int>> elements_old;
1613 for (
int ie = 0; ie < nElements_global; ie++)
1618 assert(0 <= ne && ne < nElements_global && elementFile.good());
1619 for (
int iv = 0; iv < simplexDim; iv++)
1623 assert(0 <= nv && nv < nNodes_global);
1624 element_nodes_old[iv] = nv;
1627 for (
int iv = 0; iv < simplexDim; iv++)
1630 int nN_star = element_nodes_old[iv];
1631 bool inSubdomain=
false;
1632 if (nN_star >= nodeOffsets_old[rank] && nN_star < nodeOffsets_old[rank+1])
1635 for (
int jv = 0; jv < simplexDim; jv++)
1639 int nN_star_subdomain = nN_star-nodeOffsets_old[rank];
1640 nodeStar[nN_star_subdomain].insert(element_nodes_old[jv]);
1645 elements_old[ie] = element_nodes_old;
1649 elements_old.clear();
1650 elementFile.close();
1652 PetscLogEventEnd(read_elements_event,0,0,0,0);
1653 int repartition_nodes_event;
1654 PetscLogEventRegister(
"Repart nodes",0,&repartition_nodes_event);
1655 PetscLogEventBegin(repartition_nodes_event,0,0,0,0);
1659 int max_nNodeNeighbors_node=0;
1660 for (
int nN=0;nN<nNodes_subdomain_old;nN++)
1661 max_nNodeNeighbors_node=
max(max_nNodeNeighbors_node,
static_cast<int>(nodeStar[nN].size()));
1663 PetscBool isInitialized;
1664 PetscInitialized(&isInitialized);
1665 PetscInt *nodeNeighborsOffsets_subdomain,*nodeNeighbors_subdomain,*weights_subdomain,*vertex_weights_subdomain;
1666 PetscReal *partition_weights;
1667 PetscMalloc(
sizeof(PetscInt)*(nNodes_subdomain_old+1),&nodeNeighborsOffsets_subdomain);
1668 PetscMalloc(
sizeof(PetscInt)*(nNodes_subdomain_old*max_nNodeNeighbors_node),&nodeNeighbors_subdomain);
1669 PetscMalloc(
sizeof(PetscInt)*(nNodes_subdomain_old*max_nNodeNeighbors_node),&weights_subdomain);
1670 PetscMalloc(
sizeof(PetscInt)*(nNodes_subdomain_old),&vertex_weights_subdomain);
1671 PetscMalloc(
sizeof(PetscReal)*(size),&partition_weights);
1672 for (
int sd=0;sd<size;sd++)
1673 partition_weights[sd] = 1.0/
double(size);
1674 nodeNeighborsOffsets_subdomain[0] = 0;
1676 for (
int nN = 0,offset=0; nN < nNodes_subdomain_old; nN++)
1678 for (
auto nN_star=nodeStar[nN].begin(); nN_star!=nodeStar[nN].end(); nN_star++,offset++)
1680 nodeNeighbors_subdomain[offset] = *nN_star;
1682 nodeNeighborsOffsets_subdomain[nN+1]=offset;
1685 valarray<int> neighbors(nodeNeighborsOffsets_subdomain[nN+1]-nodeNeighborsOffsets_subdomain[nN]);
1686 for (
int i=0;i<neighbors.size();i++)
1687 neighbors[i] = nodeNeighbors_subdomain[nodeNeighborsOffsets_subdomain[nN]+i];
1688 sort(&nodeNeighbors_subdomain[nodeNeighborsOffsets_subdomain[nN]],&nodeNeighbors_subdomain[nodeNeighborsOffsets_subdomain[nN+1]]);
1689 for (
int i=0;i<neighbors.size();i++)
1690 assert(neighbors[i] == nodeNeighbors_subdomain[nodeNeighborsOffsets_subdomain[nN]+i]);
1693 int weight= (nodeNeighborsOffsets_subdomain[nN+1] - nodeNeighborsOffsets_subdomain[nN]);
1694 vertex_weights_subdomain[nN] = weight;
1695 for (
int k=nodeNeighborsOffsets_subdomain[nN];k<nodeNeighborsOffsets_subdomain[nN+1];k++)
1696 weights_subdomain[k] = weight;
1702 logEvent(
"Constructing new nodal partition",5);
1704 ierr = MatCreateMPIAdj(PROTEUS_COMM_WORLD,
1705 nNodes_subdomain_old,
1707 nodeNeighborsOffsets_subdomain,
1708 nodeNeighbors_subdomain,
1710 &petscAdjacency);CHKERRABORT(PROTEUS_COMM_WORLD, ierr);
1711 const double max_rss_gb(memHardLimit/1024.0);
1712 ierr =
enforceMemoryLimit(PROTEUS_COMM_WORLD, rank, max_rss_gb,
"Done allocating MPIAdj");CHKERRABORT(PROTEUS_COMM_WORLD, ierr);
1713 MatPartitioning petscPartition;
1714 ierr = MatPartitioningCreate(PROTEUS_COMM_WORLD,&petscPartition);CHKERRABORT(PROTEUS_COMM_WORLD, ierr);
1715 ierr = MatPartitioningSetAdjacency(petscPartition,petscAdjacency);CHKERRABORT(PROTEUS_COMM_WORLD, ierr);
1716 ierr = MatPartitioningSetFromOptions(petscPartition);CHKERRABORT(PROTEUS_COMM_WORLD, ierr);
1717 ierr = MatPartitioningSetVertexWeights(petscPartition,vertex_weights_subdomain);CHKERRABORT(PROTEUS_COMM_WORLD, ierr);
1718 ierr = MatPartitioningSetPartitionWeights(petscPartition,partition_weights);CHKERRABORT(PROTEUS_COMM_WORLD, ierr);
1720 IS nodePartitioningIS_new;
1721 ierr = MatPartitioningApply(petscPartition,&nodePartitioningIS_new);CHKERRABORT(PROTEUS_COMM_WORLD, ierr);
1722 ierr = MatPartitioningDestroy(&petscPartition);CHKERRABORT(PROTEUS_COMM_WORLD, ierr);
1723 ierr =
enforceMemoryLimit(PROTEUS_COMM_WORLD, rank, max_rss_gb,
"Done applying partition");CHKERRABORT(PROTEUS_COMM_WORLD, ierr);
1726 valarray<int> nNodes_subdomain_new(size);
1727 ISPartitioningCount(nodePartitioningIS_new,size,&nNodes_subdomain_new[0]);
1729 valarray<int> nodeOffsets_new(size+1);
1730 nodeOffsets_new[0] = 0;
1731 nNodes_subdomain_max = 0;
1732 for (
int sdN = 0; sdN < size; sdN++)
1734 nodeOffsets_new[sdN+1] = nodeOffsets_new[sdN] + nNodes_subdomain_new[sdN];
1735 nNodes_subdomain_max=
max(nNodes_subdomain_max, nNodes_subdomain_new[sdN]);
1737 sprintf(log_buffer,
"Final partitioning: Max nNodes_subdomain %d nNodes_global %d",nNodes_subdomain_max,nNodes_global);
1738 logEvent(log_buffer,5);
1740 IS nodeNumberingIS_subdomain_old2new;
1741 ISPartitioningToNumbering(nodePartitioningIS_new,&nodeNumberingIS_subdomain_old2new);
1742 ierr=ISDestroy(&nodePartitioningIS_new);CHKERRABORT(PROTEUS_COMM_WORLD, ierr);
1743 const PetscInt* nodes_subdomain_old2new_array;
1744 ISGetIndices(nodeNumberingIS_subdomain_old2new, &nodes_subdomain_old2new_array);
1749 logEvent(
"Writing/reading node numberings to hdf5",5);
1750 hid_t mappings_plist_id = H5Pcreate(H5P_FILE_ACCESS);
1751 H5Pset_fapl_mpio(mappings_plist_id, PROTEUS_COMM_WORLD, MPI_INFO_NULL);
1752 const char* H5FILE_NAME(
"mappings.h5");
1753 hid_t file_id = H5Fcreate(H5FILE_NAME, H5F_ACC_TRUNC, H5P_DEFAULT, mappings_plist_id);
1754 assert(file_id != H5I_INVALID_HID);
1755 H5Pclose(mappings_plist_id);
1759 hsize_t nodes_global_count[]={
static_cast<hsize_t
>(nNodes_global)};
1760 const hsize_t ARRAY_RANK(1);
1761 hid_t nodes_old2new_filespace_id = H5Screate_simple(ARRAY_RANK, nodes_global_count, NULL);
1762 hid_t nodes_old2new_dataset_id = H5Dcreate(file_id,
"nodeNumbering_old2new", H5T_NATIVE_INT,
1763 nodes_old2new_filespace_id,
1764 H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
1765 hsize_t nodes_subdomain_count[]={
static_cast<hsize_t
>(nNodes_subdomain_old)};
1766 hsize_t nodes_subdomain_offset[]={
static_cast<hsize_t
>(nodeOffsets_old[rank])};
1767 hid_t nodes_old2new_memspace_id = H5Screate_simple(ARRAY_RANK, nodes_subdomain_count, NULL);
1768 H5Sselect_hyperslab(nodes_old2new_filespace_id, H5S_SELECT_SET,
1769 nodes_subdomain_offset, NULL,
1770 nodes_subdomain_count, NULL);
1771 hid_t nodes_old2new_plist_id = H5Pcreate(H5P_DATASET_XFER);
1772 herr_t status = H5Pset_dxpl_mpio(nodes_old2new_plist_id, H5FD_MPIO_COLLECTIVE);
1773 status = H5Dwrite(nodes_old2new_dataset_id, H5T_NATIVE_INT,
1774 nodes_old2new_memspace_id, nodes_old2new_filespace_id,
1775 nodes_old2new_plist_id, nodes_subdomain_old2new_array);
1776 H5Pclose(nodes_old2new_plist_id);
1777 H5Dclose(nodes_old2new_dataset_id);
1778 H5Sclose(nodes_old2new_memspace_id);
1779 H5Sclose(nodes_old2new_filespace_id);
1783 hid_t nodes_new2old_filespace_id = H5Screate_simple(ARRAY_RANK, nodes_global_count, NULL);
1784 hid_t nodes_new2old_dataset_id = H5Dcreate(file_id,
"nodeNumbering_new2old", H5T_NATIVE_INT,
1785 nodes_new2old_filespace_id,
1786 H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
1787 hid_t nodes_new2old_memspace_id = H5Screate_simple(ARRAY_RANK, nodes_subdomain_count, NULL);
1790 valarray<hsize_t> new_node_indices(nNodes_subdomain_old);
1791 valarray<int> old_node_indices(nNodes_subdomain_old);
1792 for (
int i=0;i<nNodes_subdomain_old;i++)
1794 new_node_indices[i] =
static_cast<hsize_t
>(nodes_subdomain_old2new_array[i]);
1795 old_node_indices[i] = nodeOffsets_old[rank]+i;
1797 status = H5Sselect_elements(nodes_new2old_filespace_id, H5S_SELECT_SET,
1798 nNodes_subdomain_old, &new_node_indices[0]);
1799 hid_t nodes_new2old_plist_id = H5Pcreate(H5P_DATASET_XFER);
1800 status = H5Pset_dxpl_mpio(nodes_new2old_plist_id, H5FD_MPIO_COLLECTIVE);
1801 status = H5Dwrite(nodes_new2old_dataset_id, H5T_NATIVE_INT,
1802 nodes_new2old_memspace_id, nodes_new2old_filespace_id,
1803 nodes_new2old_plist_id, &old_node_indices[0]);
1804 H5Pclose(nodes_new2old_plist_id);
1805 H5Dclose(nodes_new2old_dataset_id);
1806 H5Sclose(nodes_new2old_memspace_id);
1807 H5Sclose(nodes_new2old_filespace_id);
1812 IS nodeNumberingIS_global_old2new;
1813 const PetscInt *nodeNumbering_global_old2new;
1814 valarray<int> nodeNumbering_global_new2old;
1817 ISAllGather(nodeNumberingIS_subdomain_old2new,&nodeNumberingIS_global_old2new);
1818 ISGetIndices(nodeNumberingIS_global_old2new,&nodeNumbering_global_old2new);
1819 nodeNumbering_global_new2old.resize(nNodes_global);
1820 for (
int nN = 0; nN < nNodes_global; nN++)
1821 nodeNumbering_global_new2old[nodeNumbering_global_old2new[nN]] = nN;
1823 valarray<int> nodeNumbering_old2new_read(nNodes_global);
1824 file_id = H5Fopen(H5FILE_NAME, H5F_ACC_RDONLY, H5P_DEFAULT);
1825 assert(file_id != H5I_INVALID_HID);
1826 hid_t nodeNumbering_old2new_dataset_id = H5Dopen1(file_id,
"/nodeNumbering_old2new");
1827 status = H5Dread(nodeNumbering_old2new_dataset_id, H5T_NATIVE_INT, H5S_ALL, H5S_ALL, H5P_DEFAULT,
1828 &nodeNumbering_old2new_read[0]);
1829 status = H5Dclose(nodeNumbering_old2new_dataset_id);
1830 for (
int i=0;i<nNodes_global;i++)
1832 assert(nodeNumbering_global_old2new[i] == nodeNumbering_old2new_read[i]);
1834 std::cout<<
"==================out of core old2new nodes is correct!===================="<<std::endl;
1835 hid_t nodeNumbering_new2old_dataset_id = H5Dopen2(file_id,
"/nodeNumbering_new2old", H5P_DEFAULT);
1836 valarray<int> nodeNumbering_new2old_read(nNodes_global);
1837 status = H5Dread(nodeNumbering_new2old_dataset_id, H5T_NATIVE_INT, H5S_ALL, H5S_ALL, H5P_DEFAULT,
1838 &nodeNumbering_new2old_read[0]);
1839 status = H5Dclose(nodeNumbering_new2old_dataset_id);
1841 for (
int i=0;i<nNodes_global;i++)
1843 assert(nodeNumbering_global_new2old[i] == nodeNumbering_new2old_read[i]);
1845 std::cout<<
"==================out of core new2old nodes is correct!===================="<<std::endl;
1846 nodeNumbering_global_new2old.resize(0);
1849 Vec nodeNumbering_old2new_petsc;
1850 ierr = VecCreateMPI(PROTEUS_COMM_WORLD, nNodes_subdomain_old, nNodes_global, &nodeNumbering_old2new_petsc);
1851 CHKERRABORT(PROTEUS_COMM_WORLD, ierr);
1852 valarray<PetscInt> indices(nNodes_subdomain_old);
1853 valarray<PetscInt> ix(nNodes_subdomain_old);
1854 valarray<PetscScalar> y(nNodes_subdomain_old);
1855 for (
int i=0;i<nNodes_subdomain_old;i++)
1857 indices[i] = i+nodeOffsets_old[rank];
1859 y[i]=
static_cast<PetscScalar
>(nodes_subdomain_old2new_array[i]);
1861 ISLocalToGlobalMapping mapping;
1863 ierr = ISLocalToGlobalMappingCreate(PROTEUS_COMM_WORLD, bs, nNodes_subdomain_old, &indices[0], PETSC_USE_POINTER, &mapping);
1864 CHKERRABORT(PROTEUS_COMM_WORLD, ierr);
1865 ierr = VecSetLocalToGlobalMapping(nodeNumbering_old2new_petsc, mapping);
1866 CHKERRABORT(PROTEUS_COMM_WORLD, ierr);
1867 ierr = VecSetValuesLocal(nodeNumbering_old2new_petsc, nNodes_subdomain_old, &ix[0], &y[0], INSERT_VALUES);
1868 CHKERRABORT(PROTEUS_COMM_WORLD, ierr);
1869 ierr = VecAssemblyBegin(nodeNumbering_old2new_petsc);
1870 CHKERRABORT(PROTEUS_COMM_WORLD, ierr);
1871 ierr = VecAssemblyEnd(nodeNumbering_old2new_petsc);
1872 CHKERRABORT(PROTEUS_COMM_WORLD, ierr);
1875 PetscScalar* test_local;
1876 VecGetArray(nodeNumbering_old2new_petsc,&test_local);
1877 for (
int i=0;i<nNodes_subdomain_old;i++)
1879 assert(test_local[i] ==
static_cast<PetscScalar
>(nodes_subdomain_old2new_array[i]));
1881 VecRestoreArray(nodeNumbering_old2new_petsc,&test_local);
1886 PetscLogEventEnd(repartition_nodes_event,0,0,0,0);
1887 int receive_element_mask_event;
1888 PetscLogEventRegister(
"Recv. ele mask",0,&receive_element_mask_event);
1889 PetscLogEventBegin(receive_element_mask_event,0,0,0,0);
1894 logEvent(
"Collecting elements containing locally owned nodes (element support of the subdomain nodes)",5);
1895 PetscLogEventEnd(receive_element_mask_event,0,0,0,0);
1896 ierr =
enforceMemoryLimit(PROTEUS_COMM_WORLD, rank, max_rss_gb,
"Done with masks");CHKERRABORT(PROTEUS_COMM_WORLD, ierr);
1897 int build_subdomains_reread_elements_event;
1898 PetscLogEventRegister(
"Reread eles",0,&build_subdomains_reread_elements_event);
1899 PetscLogEventBegin(build_subdomains_reread_elements_event,0,0,0,0);
1906 std::ifstream elementFile2(elementFileName.c_str());
1907 if (!elementFile2.good())
1909 std::cerr<<
"cannot open Tetgen elements file"
1910 <<elementFileName<<std::endl;
1914 elementFile2 >>
eatcomments >> nElements_global >> nNodesPerSimplex >> hasElementMarkers >>
eatline;
1915 assert(nElements_global > 0);
1916 assert(nNodesPerSimplex == simplexDim);
1918 set<int> elements_subdomain_owned;
1919 valarray<int> element_nodes_new(4);
1922 map<int,valarray<int>> elementNodesArrayMap;
1923 map<int,long int> elementMaterialTypesMap;
1925 map<NodeTuple<2>,set<pair<int,int>>> edgeElementsMap;
1926 set<int> node_collection;
1930 map<int, int> nodes_old2new_subdomain_map;
1933 for (
int ie = 0; ie < nElements_global; ie++)
1942 int ne, nv, elementId(0);
1943 long double elementId_double;
1946 assert(0 <= ne && ne < nElements_global && elementFile.good());
1947 elements_collection.push_back(ne);
1948 for (
int iv = 0; iv < simplexDim; iv++)
1950 elementFile2 >> nv ;
1952 assert(0 <= nv && nv < nNodes_global);
1953 node_collection.insert(nv);
1954 element_nodes_old[iv] = nv;
1956 if (hasElementMarkers > 0)
1958 elementFile2 >> elementId_double;
1959 elementId_collection.push_back(elementId_double);
1961 element_old_nodes_collection.push_back(element_nodes_old);
1962 if (elements_collection.size() == nElements_global/size || ie == nElements_global-1)
1965 valarray<int> nodes_old2new_subset(node_collection.size());
1966 map<int,int> nodes_old2new_subset_map;
1967 bool test_scatter=
false;
1970 Vec node_old2new_collection_petsc;
1971 ierr = VecCreateMPI(PROTEUS_COMM_WORLD, node_collection.size(), PETSC_DETERMINE,&node_old2new_collection_petsc);
1972 CHKERRABORT(PROTEUS_COMM_WORLD, ierr);
1973 PetscInt collection_low, collection_high;
1974 ierr = VecGetOwnershipRange(node_old2new_collection_petsc, &collection_low, &collection_high);
1975 CHKERRABORT(PROTEUS_COMM_WORLD, ierr);
1976 valarray<PetscInt> node_collection_array(node_collection.size());
1977 valarray<PetscInt> iy_collection(node_collection.size());
1979 for (
auto nv=node_collection.begin();nv!=node_collection.end();nv++,i_nc++)
1981 node_collection_array[i_nc] =
static_cast<PetscInt
>(*nv);
1982 iy_collection[i_nc] = i_nc + collection_low;
1985 ierr = ISCreateGeneral(PROTEUS_COMM_WORLD, node_collection.size(), &node_collection_array[0],
1986 PETSC_USE_POINTER, &old_ix);
1987 CHKERRABORT(PROTEUS_COMM_WORLD, ierr);
1988 ierr = ISCreateGeneral(PROTEUS_COMM_WORLD, node_collection.size(), &iy_collection[0],
1989 PETSC_USE_POINTER, &iy);
1990 CHKERRABORT(PROTEUS_COMM_WORLD, ierr);
1991 VecScatter old2new_collection_scatter;
1992 ierr = VecScatterCreate(nodeNumbering_old2new_petsc, old_ix,
1993 node_old2new_collection_petsc, iy,
1994 &old2new_collection_scatter);
1995 CHKERRABORT(PROTEUS_COMM_WORLD, ierr);
1996 ierr = VecScatterBegin(old2new_collection_scatter,
1997 nodeNumbering_old2new_petsc, node_old2new_collection_petsc,
1998 INSERT_VALUES, SCATTER_FORWARD);
1999 CHKERRABORT(PROTEUS_COMM_WORLD, ierr);
2001 ierr = VecScatterEnd(old2new_collection_scatter,
2002 nodeNumbering_old2new_petsc, node_old2new_collection_petsc,
2003 INSERT_VALUES, SCATTER_FORWARD);
2004 CHKERRABORT(PROTEUS_COMM_WORLD, ierr);
2005 PetscScalar* nodes_old2new_subset_array;
2006 ierr = VecGetArray(node_old2new_collection_petsc, &nodes_old2new_subset_array);
2007 CHKERRABORT(PROTEUS_COMM_WORLD, ierr);
2008 for (
int i=0;i<node_collection.size();i++)
2010 nodes_old2new_subset[i] = nodes_old2new_subset_array[i];
2012 ierr = VecRestoreArray(node_old2new_collection_petsc, &nodes_old2new_subset_array);
2013 CHKERRABORT(PROTEUS_COMM_WORLD, ierr);
2014 ierr = VecDestroy(&node_old2new_collection_petsc);
2015 CHKERRABORT(PROTEUS_COMM_WORLD, ierr);
2016 ierr = VecScatterDestroy(&old2new_collection_scatter);
2017 CHKERRABORT(PROTEUS_COMM_WORLD, ierr);
2018 ierr = ISDestroy(&old_ix);
2019 CHKERRABORT(PROTEUS_COMM_WORLD, ierr);
2020 ierr = ISDestroy(&iy);
2021 CHKERRABORT(PROTEUS_COMM_WORLD, ierr);
2022 for (
int i=0;i<node_collection.size();i++)
2024 nodes_old2new_subset_map[node_collection_array[i]] = nodes_old2new_subset[i];
2028 for (
int i=0;i<node_collection.size();i++)
2030 assert(nodes_old2new_subset_map[node_collection_array[i]] == nodeNumbering_global_old2new[node_collection_array[i]]);
2037 valarray<hsize_t> node_collection_array(node_collection.size());
2039 for (
auto nv=node_collection.begin();nv!=node_collection.end();nv++,i_nc++)
2041 node_collection_array[i_nc] =
static_cast<hsize_t
>(*nv);
2043 file_id = H5Fopen(H5FILE_NAME, H5F_ACC_RDONLY, H5P_DEFAULT);
2044 assert(file_id != H5I_INVALID_HID);
2045 nodes_old2new_dataset_id = H5Dopen1(file_id,
"/nodeNumbering_old2new");
2046 nodes_old2new_filespace_id = H5Dget_space(nodes_old2new_dataset_id);
2047 status = H5Sselect_elements(nodes_old2new_filespace_id, H5S_SELECT_SET,
2048 node_collection.size(), &node_collection_array[0]);
2049 hsize_t dims[] = {
static_cast<hsize_t
>(node_collection.size())};
2050 hid_t nodes_old2new_subset_memspace_id = H5Screate_simple(1, dims, NULL);
2051 nodes_old2new_plist_id = H5Pcreate(H5P_DATASET_XFER);
2052 status = H5Pset_dxpl_mpio(nodes_old2new_plist_id, H5FD_MPIO_COLLECTIVE);
2053 status = H5Dread(nodes_old2new_dataset_id, H5T_NATIVE_INT,
2054 nodes_old2new_subset_memspace_id, nodes_old2new_filespace_id,
2055 H5P_DEFAULT, &nodes_old2new_subset[0]);
2056 H5Pclose(nodes_old2new_plist_id);
2057 H5Sclose(nodes_old2new_subset_memspace_id);
2058 H5Sclose(nodes_old2new_filespace_id);
2059 H5Dclose(nodes_old2new_dataset_id);
2061 for (
int i=0;i<node_collection.size();i++)
2063 nodes_old2new_subset_map[node_collection_array[i]] = nodes_old2new_subset[i];
2067 for (
int i=0;i<node_collection.size();i++)
2069 assert(nodes_old2new_subset_map[node_collection_array[i]] == nodeNumbering_global_old2new[node_collection_array[i]]);
2073 for (
int eN_c = eN_c_start; eN_c < ie+1;eN_c++)
2075 ne = elements_collection[eN_c-eN_c_start];
2076 for (
int iv = 0; iv < simplexDim; iv++)
2078 element_nodes_old[iv] = element_old_nodes_collection[eN_c-eN_c_start][iv];
2079 element_nodes_new[iv] = nodes_old2new_subset_map[element_nodes_old[iv]];
2082 for (
int iv = 0; iv < simplexDim; iv++)
2084 int nN_star_new = element_nodes_new[iv];
2085 bool inSubdomain=
false;
2086 if (nN_star_new >= nodeOffsets_new[rank] && nN_star_new < nodeOffsets_new[rank+1])
2089 for (
int iv=0; iv < simplexDim; iv++)
2090 nodes_old2new_subdomain_map[element_nodes_old[iv]] = element_nodes_new[iv];
2092 for (
int ebN=0;ebN < 4 ; ebN++)
2094 int nodes[3] = { element_nodes_new[(ebN+1) % 4],
2095 element_nodes_new[(ebN+2) % 4],
2096 element_nodes_new[(ebN+3) % 4]};
2098 if(elementBoundaryElementsMap.find(nodeTuple) != elementBoundaryElementsMap.end())
2100 if (elementBoundaryElementsMap[nodeTuple].right == -1 && ne != elementBoundaryElementsMap[nodeTuple].left)
2102 elementBoundaryElementsMap[nodeTuple].right=ne;
2103 elementBoundaryElementsMap[nodeTuple].right_ebN_element=ebN;
2112 for (
int nNL=0,edN=0;nNL < 4 ; nNL++)
2113 for(
int nNR=nNL+1;nNR < 4;nNR++,edN++)
2115 int nodes[2] = { element_nodes_new[nNL],
2116 element_nodes_new[nNR]};
2118 edgeElementsMap[nodeTuple].insert(pair<int,int>(ne,edN));
2121 int nN_star_new_subdomain = nN_star_new - nodeOffsets_new[rank];
2122 nodeElementsStar[nN_star_new_subdomain].insert(ne);
2123 for (
int jv = 0; jv < simplexDim; jv++)
2127 int nN_point_new = element_nodes_new[jv];
2128 nodeStarNew[nN_star_new_subdomain].insert(nN_point_new);
2134 elementNodesArrayMap[ne] = element_nodes_new;
2138 if (elementNodesArrayMap.find(ne) != elementNodesArrayMap.end())
2140 if (nodeTuple.
nodes[1] >= nodeOffsets_new[rank] && nodeTuple.
nodes[1] < nodeOffsets_new[rank+1])
2141 elements_subdomain_owned.insert(ne);
2142 if (hasElementMarkers > 0)
2144 elementId =
static_cast<long int>(elementId_collection[eN_c-eN_c_start]);
2145 elementMaterialTypesMap[ne] = elementId;
2150 node_collection.clear();
2151 elements_collection.clear();
2152 elementId_collection.clear();
2153 element_old_nodes_collection.clear();
2157 elementFile2.close();
2158 int nElements_owned_subdomain(elements_subdomain_owned.size()),
2159 nElements_owned_new=0;
2160 MPI_Allreduce(&nElements_owned_subdomain,&nElements_owned_new,1,MPI_INT,MPI_SUM,PROTEUS_COMM_WORLD);
2161 assert(nElements_owned_new == nElements_global);
2162 PetscLogEventEnd(build_subdomains_reread_elements_event,0,0,0,0);
2163 int build_subdomains_send_marked_elements_event;
2164 PetscLogEventRegister(
"Mark/send eles",0,&build_subdomains_send_marked_elements_event);
2165 PetscLogEventBegin(build_subdomains_send_marked_elements_event,0,0,0,0);
2166 ierr =
enforceMemoryLimit(PROTEUS_COMM_WORLD, rank, max_rss_gb,
"Done marking elements");CHKERRABORT(PROTEUS_COMM_WORLD, ierr);
2172 valarray<int> nodeElementOffsets(nNodes_subdomain_new[rank]+1);
2173 nodeElementOffsets[0] = 0;
2174 for (
int nN = 0; nN < nNodes_subdomain_new[rank]; nN++)
2175 nodeElementOffsets[nN+1] = nodeElementOffsets[nN]+nodeElementsStar[nN].size();
2176 valarray<int> nodeElementsArray(nodeElementOffsets[nNodes_subdomain_new[rank]]);
2177 for (
int nN=0,offset=0; nN < nNodes_subdomain_new[rank]; nN++)
2179 for (set<int>::iterator eN_star = nodeElementsStar[nN].begin(); eN_star != nodeElementsStar[nN].end();
2182 nodeElementsArray[offset] = *eN_star;
2185 nodeElementsStar.clear();
2187 valarray<int> nodeStarOffsetsNew(nNodes_subdomain_new[rank]+1);
2188 nodeStarOffsetsNew[0] = 0;
2189 for (
int nN=1;nN<nNodes_subdomain_new[rank]+1;nN++)
2190 nodeStarOffsetsNew[nN] = nodeStarOffsetsNew[nN-1] + nodeStarNew[nN-1].size();
2191 valarray<int> nodeStarArrayNew(nodeStarOffsetsNew[nNodes_subdomain_new[rank]]);
2192 for (
int nN=0,offset=0;nN<nNodes_subdomain_new[rank];nN++)
2194 for (set<int>::iterator nN_star=nodeStarNew[nN].begin();nN_star!=nodeStarNew[nN].end();nN_star++,offset++)
2196 nodeStarArrayNew[offset] = *nN_star;
2199 nodeStarNew.clear();
2200 PetscLogEventEnd(build_subdomains_send_marked_elements_event,0,0,0,0);
2201 int build_subdomains_global_numbering_elements_event;
2202 PetscLogEventRegister(
"Global ele nmbr",0,&build_subdomains_global_numbering_elements_event);
2203 PetscLogEventBegin(build_subdomains_global_numbering_elements_event,0,0,0,0);
2207 logEvent(
"Generating global element numbering corresponding to new subdomain ownership",5);
2208 valarray<int> nElements_subdomain_new(size),
2209 elementOffsets_new(size+1);
2210 for (
int sdN = 0; sdN < size; sdN++)
2213 nElements_subdomain_new[sdN] = int(elements_subdomain_owned.size());
2215 nElements_subdomain_new[sdN] = 0;
2217 valarray<int> nElements_subdomain_new_send = nElements_subdomain_new;
2218 MPI_Allreduce(&nElements_subdomain_new_send[0],&nElements_subdomain_new[0],size,MPI_INT,MPI_SUM,PROTEUS_COMM_WORLD);
2219 sprintf(log_buffer,
"Final partitioning: Max nElements_subdomain %d nElements_global %d",nElements_subdomain_new.max(),nElements_global);
2220 logEvent(log_buffer,5);
2223 elementOffsets_new[0] = 0;
2224 for (
int sdN = 0; sdN < size; sdN++)
2225 elementOffsets_new[sdN+1] = elementOffsets_new[sdN] + nElements_subdomain_new[sdN];
2227 valarray<int> elementNumbering_subdomain_new2old(elements_subdomain_owned.size());
2228 set<int>::iterator eN_ownedp = elements_subdomain_owned.begin();
2229 map<int, int> elementNumbering_old2new_map;
2230 for (
int eN = 0; eN < int(elements_subdomain_owned.size()); eN++,eN_ownedp++)
2233 elementNumbering_subdomain_new2old[eN] = *eN_ownedp;
2234 elementNumbering_old2new_map[*eN_ownedp] = eN;
2239 ierr =
enforceMemoryLimit(PROTEUS_COMM_WORLD, rank, max_rss_gb,
"Writing element numberings to hdf5");CHKERRABORT(PROTEUS_COMM_WORLD, ierr);
2240 logEvent(
"Writing element numberings to hdf5",5);
2241 mappings_plist_id = H5Pcreate(H5P_FILE_ACCESS);
2242 status = H5Pset_fapl_mpio(mappings_plist_id, PROTEUS_COMM_WORLD, MPI_INFO_NULL);
2243 assert(status >= 0);
2244 file_id = H5Fopen(H5FILE_NAME, H5F_ACC_RDWR, mappings_plist_id);
2245 assert(file_id != H5I_INVALID_HID);
2246 H5Pclose(mappings_plist_id);
2247 hsize_t e_dims[]={
static_cast<hsize_t
>(nElements_global)};
2248 hid_t e_new2old_filespace_id = H5Screate_simple(ARRAY_RANK, e_dims, NULL);
2249 hid_t e_new2old_dataset_id = H5Dcreate(file_id,
"elementNumbering_new2old", H5T_NATIVE_INT,
2250 e_new2old_filespace_id,
2251 H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
2252 hsize_t e_count[]={
static_cast<hsize_t
>(nElements_subdomain_new[rank])};
2253 hsize_t e_offset[]={
static_cast<hsize_t
>(elementOffsets_new[rank])};
2254 hid_t e_new2old_memspace_id = H5Screate_simple(ARRAY_RANK, e_count, NULL);
2255 H5Sselect_hyperslab(e_new2old_filespace_id, H5S_SELECT_SET, e_offset, NULL, e_count, NULL);
2256 hid_t e_new2old_plist_id = H5Pcreate(H5P_DATASET_XFER);
2257 H5Pset_dxpl_mpio(e_new2old_plist_id, H5FD_MPIO_COLLECTIVE);
2258 status = H5Dwrite(e_new2old_dataset_id, H5T_NATIVE_INT,
2259 e_new2old_memspace_id, e_new2old_filespace_id,
2260 e_new2old_plist_id, &elementNumbering_subdomain_new2old[0]);
2261 H5Dflush(e_new2old_dataset_id);
2262 H5Pclose(e_new2old_plist_id);
2263 H5Dclose(e_new2old_dataset_id);
2264 H5Sclose(e_new2old_memspace_id);
2265 H5Sclose(e_new2old_filespace_id);
2266 valarray<hsize_t> old_element_indices(
static_cast<hsize_t
>(nElements_subdomain_new[rank]));
2267 valarray<int> new_element_indices(nElements_subdomain_new[rank]);
2268 ierr =
enforceMemoryLimit(PROTEUS_COMM_WORLD, rank, max_rss_gb,
"new2old done; now old2new");CHKERRABORT(PROTEUS_COMM_WORLD, ierr);
2270 ierr =
enforceMemoryLimit(PROTEUS_COMM_WORLD, rank, max_rss_gb,
"old2new index valarrays done");CHKERRABORT(PROTEUS_COMM_WORLD, ierr);
2271 for (
int i=0;i<nElements_subdomain_new[rank];i++)
2273 old_element_indices[i] =
static_cast<hsize_t
>(elementNumbering_subdomain_new2old[i]);
2274 new_element_indices[i] = elementOffsets_new[rank]+i;
2276 hid_t e_old2new_filespace_id = H5Screate_simple(ARRAY_RANK, e_dims, NULL);
2277 hid_t e_old2new_dataset_id = H5Dcreate(file_id,
"elementNumbering_old2new", H5T_NATIVE_INT,
2278 e_old2new_filespace_id,
2279 H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
2280 hid_t e_old2new_memspace_id = H5Screate_simple(ARRAY_RANK, e_count, NULL);
2281 status = H5Sselect_elements(e_old2new_filespace_id, H5S_SELECT_SET,
2282 nElements_subdomain_new[rank], &old_element_indices[0]);
2283 hid_t e_old2new_plist_id = H5Pcreate(H5P_DATASET_XFER);
2284 H5Pset_dxpl_mpio(e_old2new_plist_id, H5FD_MPIO_COLLECTIVE);
2285 status = H5Dwrite(e_old2new_dataset_id, H5T_NATIVE_INT,
2286 e_old2new_memspace_id, e_old2new_filespace_id,
2287 e_old2new_plist_id, &new_element_indices[0]);
2288 H5Pclose(e_old2new_plist_id);
2289 H5Sclose(e_old2new_memspace_id);
2290 H5Dclose(e_old2new_dataset_id);
2291 H5Sclose(e_old2new_filespace_id);
2296 ierr =
enforceMemoryLimit(PROTEUS_COMM_WORLD, rank, max_rss_gb,
"Reading numberings from hdf5 to get subdomain map");CHKERRABORT(PROTEUS_COMM_WORLD, ierr);
2297 logEvent(
"Reading element numberings from hdf5 to get subdomain map",5);
2299 valarray<hsize_t> old_element_indices_subdomain(
static_cast<hsize_t
>(elementNodesArrayMap.size()));
2300 valarray<int> new_element_indices_subdomain(elementNodesArrayMap.size());
2301 int eN_old_subdomain = 0;
2302 for (
auto it = elementNodesArrayMap.begin(); it != elementNodesArrayMap.end(); it++, eN_old_subdomain++)
2304 old_element_indices_subdomain[eN_old_subdomain] = it->first;
2306 file_id = H5Fopen(H5FILE_NAME, H5F_ACC_RDONLY, H5P_DEFAULT);
2307 e_old2new_dataset_id = H5Dopen1(file_id,
"/elementNumbering_old2new");
2308 e_old2new_filespace_id = H5Dget_space(e_old2new_dataset_id);
2309 status = H5Sselect_elements(e_old2new_filespace_id, H5S_SELECT_SET,
2310 elementNodesArrayMap.size(), &old_element_indices_subdomain[0]);
2311 hsize_t e_subdomain_count[]={
static_cast<hsize_t
>(elementNodesArrayMap.size())};
2312 hid_t e_old2new_subdomain_memspace_id = H5Screate_simple(ARRAY_RANK, e_subdomain_count, NULL);
2313 hid_t e_old2new_subdomain_plist_id = H5Pcreate(H5P_DATASET_XFER);
2314 status = H5Pset_dxpl_mpio(e_old2new_subdomain_plist_id, H5FD_MPIO_COLLECTIVE);
2315 status = H5Dread(e_old2new_dataset_id, H5T_NATIVE_INT,
2316 e_old2new_subdomain_memspace_id, e_old2new_filespace_id,
2317 H5P_DEFAULT, &new_element_indices_subdomain[0]);
2318 H5Pclose(e_old2new_subdomain_plist_id);
2319 H5Sclose(e_old2new_subdomain_memspace_id);
2320 H5Dclose(e_old2new_dataset_id);
2321 H5Sclose(e_old2new_filespace_id);
2323 ierr =
enforceMemoryLimit(PROTEUS_COMM_WORLD, rank, max_rss_gb,
"Constructing element subdomain maps");CHKERRABORT(PROTEUS_COMM_WORLD, ierr);
2324 map<int,int> elementNumbering_old2new_subdomain_map;
2325 map<int,int> elementNumbering_new2old_subdomain_map;
2326 for (
int i=0;i<elementNodesArrayMap.size();i++)
2328 elementNumbering_old2new_subdomain_map[old_element_indices_subdomain[i]] = new_element_indices_subdomain[i];
2329 elementNumbering_new2old_subdomain_map[new_element_indices_subdomain[i]] = old_element_indices_subdomain[i];
2334 ISRestoreIndices(nodeNumberingIS_global_old2new,&nodeNumbering_global_old2new);
2335 ierr=ISDestroy(&nodeNumberingIS_global_old2new);CHKERRABORT(PROTEUS_COMM_WORLD, ierr);
2337 IS elementNumberingIS_subdomain_new2old;
2338 const PetscInt *elementNumbering_global_new2old;
2339 valarray<int> elementNumbering_global_old2new(nElements_global);
2340 ISCreateGeneral(PROTEUS_COMM_WORLD,elements_subdomain_owned.size(),
2341 &elementNumbering_subdomain_new2old[0],PETSC_COPY_VALUES,
2342 &elementNumberingIS_subdomain_new2old);
2343 IS elementNumberingIS_global_new2old;
2344 ISAllGather(elementNumberingIS_subdomain_new2old,&elementNumberingIS_global_new2old);
2345 ISGetIndices(elementNumberingIS_global_new2old,&elementNumbering_global_new2old);
2347 for (
int eN = 0; eN < nElements_global; eN++)
2349 elementNumbering_global_old2new[elementNumbering_global_new2old[eN]] = eN;
2352 valarray<int> elementNumbering_global_old2new_read(nElements_global);
2353 file_id = H5Fopen(H5FILE_NAME, H5F_ACC_RDONLY, H5P_DEFAULT);
2354 hid_t e_old2new_dataset_id = H5Dopen1(file_id,
"/elementNumbering_old2new");
2355 status = H5Dread(e_old2new_dataset_id, H5T_NATIVE_INT, H5S_ALL, H5S_ALL, H5P_DEFAULT,
2356 &elementNumbering_global_old2new_read[0]);
2357 status = H5Dclose(e_old2new_dataset_id);
2358 for (
int i=0;i<nElements_global;i++)
2360 assert(elementNumbering_global_old2new[i] == elementNumbering_global_old2new_read[i]);
2362 std::cout<<
"==================out of core elements old2new is correct!===================="<<std::endl;
2363 for (
auto it = elementNodesArrayMap.begin(); it != elementNodesArrayMap.end(); it++)
2364 assert(elementNumbering_old2new_subdomain_map[it->first] ==
2365 elementNumbering_global_old2new_read[it->first]);
2366 std::cout<<
"==================out of core elements old2new map is correct!===================="<<std::endl;
2368 valarray<int> elementNumbering_global_new2old_read(nElements_global);
2369 hid_t e_new2old_dataset_id = H5Dopen1(file_id,
"/elementNumbering_new2old");
2370 status = H5Dread(e_new2old_dataset_id, H5T_NATIVE_INT, H5S_ALL, H5S_ALL, H5P_DEFAULT,
2371 &elementNumbering_global_new2old_read[0]);
2372 status = H5Dclose(e_new2old_dataset_id);
2373 status = H5Fclose(file_id);
2374 for (
int i=0;i<nElements_global;i++)
2375 assert(elementNumbering_global_new2old[i] == elementNumbering_global_new2old_read[i]);
2376 std::cout<<
"==================out of core elements new2old is correct!===================="<<std::endl;
2377 for (
size_t i=0;i< new_element_indices_subdomain.size(); i++)
2379 int eN = new_element_indices_subdomain[i];
2380 assert(elementNumbering_new2old_subdomain_map[eN] == elementNumbering_global_new2old_read[eN]);
2382 std::cout<<
"==================out of core elements new2old map is correct!===================="<<std::endl;
2383 ISRestoreIndices(elementNumberingIS_global_new2old,&elementNumbering_global_new2old);
2384 ierr=ISDestroy(&elementNumberingIS_global_new2old);CHKERRABORT(PROTEUS_COMM_WORLD, ierr);
2385 ierr=ISDestroy(&elementNumberingIS_subdomain_new2old);CHKERRABORT(PROTEUS_COMM_WORLD, ierr);
2387 PetscLogEventEnd(build_subdomains_global_numbering_elements_event,0,0,0,0);
2388 ierr =
enforceMemoryLimit(PROTEUS_COMM_WORLD, rank, max_rss_gb,
"Done allocating element numbering new2old/old2new");CHKERRABORT(PROTEUS_COMM_WORLD, ierr);
2389 int build_subdomains_faces_event;
2390 PetscLogEventRegister(
"Subd faces",0,&build_subdomains_faces_event);
2391 PetscLogEventBegin(build_subdomains_faces_event,0,0,0,0);
2398 logEvent(
"Generating global face numbering and ownership corresponding to new subdomain ownership",5);
2399 std::ifstream elementBoundaryFile(elementBoundaryFileName.c_str());
2401 if (!elementBoundaryFile.good())
2403 std::cerr<<
"cannot open Tetgen face file "
2404 <<elementBoundaryFileName<<std::endl;
2409 bool hasElementBoundaryMarkers =
false;
2410 int nElementBoundaries_global;
2411 int ihasElementBoundaryMarkers(0);
2412 elementBoundaryFile >>
eatcomments >> nElementBoundaries_global >> ihasElementBoundaryMarkers >>
eatline ;
2413 assert(nElementBoundaries_global > 0);
2414 if (ihasElementBoundaryMarkers > 0)
2416 hasElementBoundaryMarkers =
true;
2420 set<int> elementBoundaries_subdomain_owned;
2421 vector<set<int>> nodeElementBoundariesStar(nNodes_subdomain_new[rank]);
2422 map<int,int> elementBoundaryMaterialTypesMap;
2423 map<int,valarray<int>> elementBoundariesMap;
2424 set<int> supportedElementBoundaries;
2425 set<int> elementBoundaries_subdomain;
2426 for (
int ieb = 0; ieb < nElementBoundaries_global; ieb++)
2428 int neb,nn0,nn1,nn2;
int ebId(0);
2429 elementBoundaryFile >>
eatcomments >> neb >> nn0 >> nn1 >> nn2;
2430 if (ihasElementBoundaryMarkers > 0)
2431 elementBoundaryFile >> ebId;
2436 assert(0 <= neb && neb < nElementBoundaries_global && elementBoundaryFile.good());
2441 if (nodes_old2new_subdomain_map.find(nn0) != nodes_old2new_subdomain_map.end())
2442 nn0_new = nodes_old2new_subdomain_map[nn0];
2443 if (nn0_new >= nodeOffsets_new[rank] && nn0_new < nodeOffsets_new[rank+1])
2445 nodeElementBoundariesStar[nn0_new-nodeOffsets_new[rank]].insert(neb);
2446 supportedElementBoundaries.insert(neb);
2449 if (nodes_old2new_subdomain_map.find(nn1) != nodes_old2new_subdomain_map.end())
2450 nn1_new = nodes_old2new_subdomain_map[nn1];
2451 if (nn1_new >= nodeOffsets_new[rank] && nn1_new < nodeOffsets_new[rank+1])
2453 nodeElementBoundariesStar[nn1_new-nodeOffsets_new[rank]].insert(neb);
2454 supportedElementBoundaries.insert(neb);
2457 if (nodes_old2new_subdomain_map.find(nn2) != nodes_old2new_subdomain_map.end())
2458 nn2_new = nodes_old2new_subdomain_map[nn2];
2459 if (nn2_new >= nodeOffsets_new[rank] && nn2_new < nodeOffsets_new[rank+1])
2461 nodeElementBoundariesStar[nn2_new-nodeOffsets_new[rank]].insert(neb);
2462 supportedElementBoundaries.insert(neb);
2466 const int nodes[3] = {nn0_new,nn1_new,nn2_new};
2468 elementBoundaryFile >>
eatline;
2469 if (elementBoundaryElementsMap.find(nodeTuple) != elementBoundaryElementsMap.end())
2471 elementBoundaries_subdomain.insert(neb);
2472 assert(nn0_new >= 0 && nn1_new >= 0 && nn2_new >= 0);
2473 if (nodeTuple.
nodes[1] >= nodeOffsets_new[rank] && nodeTuple.
nodes[1] < nodeOffsets_new[rank+1])
2474 elementBoundaries_subdomain_owned.insert(neb);
2475 if (ihasElementBoundaryMarkers > 0)
2476 elementBoundaryMaterialTypesMap[neb]=ebId;
2477 int eN_left = elementNumbering_old2new_subdomain_map[elementBoundaryElementsMap[nodeTuple].left];
2478 if (elementBoundariesMap.find(eN_left) != elementBoundariesMap.end())
2480 elementBoundariesMap[eN_left][elementBoundaryElementsMap[nodeTuple].left_ebN_element] = neb;
2485 valarray<int> elementBoundaries_element(-1,4);
2486 elementBoundariesMap[eN_left] = elementBoundaries_element;
2488 elementBoundariesMap[eN_left][elementBoundaryElementsMap[nodeTuple].left_ebN_element] = neb;
2490 if (elementBoundaryElementsMap[nodeTuple].right >= 0)
2492 int eN_right = elementNumbering_old2new_subdomain_map[elementBoundaryElementsMap[nodeTuple].right];
2493 if (elementBoundariesMap.find(eN_right) != elementBoundariesMap.end())
2495 elementBoundariesMap[eN_right][elementBoundaryElementsMap[nodeTuple].right_ebN_element] = neb;
2500 valarray<int> elementBoundaries_element(-1, 4);
2501 elementBoundariesMap[eN_right] = elementBoundaries_element;
2503 elementBoundariesMap[eN_right][elementBoundaryElementsMap[nodeTuple].right_ebN_element] = neb;
2509 elementBoundaryFile.close();
2510 ierr =
enforceMemoryLimit(PROTEUS_COMM_WORLD, rank, max_rss_gb,
"Done reading element boundaries");CHKERRABORT(PROTEUS_COMM_WORLD, ierr);
2511 int nElementBoundaries_owned_subdomain=elementBoundaries_subdomain_owned.size(),
2512 nElementBoundaries_owned_new=0;
2513 MPI_Allreduce(&nElementBoundaries_owned_subdomain,&nElementBoundaries_owned_new,1,MPI_INT,MPI_SUM,PROTEUS_COMM_WORLD);
2514 assert(nElementBoundaries_owned_new == nElementBoundaries_global);
2516 for (
auto elementBoundariesp=elementBoundariesMap.begin();
2517 elementBoundariesp!=elementBoundariesMap.end();
2518 elementBoundariesp++)
2521 for (
int iv=0;iv<4;iv++)
2524 int nN_global = elementNodesArrayMap[elementNumbering_new2old_subdomain_map[elementBoundariesp->first]][iv];
2525 if (nN_global >= nodeOffsets_new[rank] && nN_global < nodeOffsets_new[rank+1])
2528 for(
int eb=0;eb<4;eb++)
2530 nodeElementBoundariesStar[nN_global-nodeOffsets_new[rank]].insert(elementBoundariesp->second[eb]);
2536 valarray<int> nodeElementBoundaryOffsets(nNodes_subdomain_new[rank]+1);
2537 nodeElementBoundaryOffsets[0] = 0;
2538 for (
int nN = 0; nN < nNodes_subdomain_new[rank]; nN++)
2539 nodeElementBoundaryOffsets[nN+1] = nodeElementBoundaryOffsets[nN]+nodeElementBoundariesStar[nN].size();
2540 valarray<int> nodeElementBoundariesArray(nodeElementBoundaryOffsets[nNodes_subdomain_new[rank]]);
2541 for (
int nN=0,offset=0; nN < nNodes_subdomain_new[rank]; nN++)
2543 for (set<int>::iterator ebN_star = nodeElementBoundariesStar[nN].begin(); ebN_star != nodeElementBoundariesStar[nN].end();
2544 ebN_star++,offset++)
2546 nodeElementBoundariesArray[offset] = *ebN_star;
2551 valarray<int> nElementBoundaries_subdomain_new(size),
2552 elementBoundaryOffsets_new(size+1);
2553 for (
int sdN=0;sdN<size;sdN++)
2555 nElementBoundaries_subdomain_new[sdN] = elementBoundaries_subdomain_owned.size();
2557 nElementBoundaries_subdomain_new[sdN] = 0;
2558 valarray<int> nElementBoundaries_subdomain_new_send=nElementBoundaries_subdomain_new;
2559 MPI_Allreduce(&nElementBoundaries_subdomain_new_send[0],&nElementBoundaries_subdomain_new[0],
2560 size,MPI_INT,MPI_SUM,PROTEUS_COMM_WORLD);
2561 elementBoundaryOffsets_new[0] = 0;
2562 for (
int sdN=0;sdN<size;sdN++)
2563 elementBoundaryOffsets_new[sdN+1] = elementBoundaryOffsets_new[sdN]+nElementBoundaries_subdomain_new[sdN];
2567 valarray<int> elementBoundaryNumbering_new2old(elementBoundaries_subdomain_owned.size());
2568 map<int,int> elementBoundaryNumbering_old2new_map;
2569 set<int>::iterator ebN_ownedp=elementBoundaries_subdomain_owned.begin();
2570 for (
int ebN=0;ebN<int(elementBoundaries_subdomain_owned.size());ebN++, ebN_ownedp++)
2572 elementBoundaryNumbering_new2old[ebN] = *ebN_ownedp;
2573 elementBoundaryNumbering_old2new_map[*ebN_ownedp] = ebN+elementBoundaryOffsets_new[rank];
2578 logEvent(
"Writing/reading elementBoundary (face) numberings to hdf5",5);
2579 hsize_t eb_dims[]={
static_cast<hsize_t
>(nElementBoundaries_global)};
2580 hid_t eb_new2old_filespace_id = H5Screate_simple(ARRAY_RANK, eb_dims, NULL);
2581 mappings_plist_id = H5Pcreate(H5P_FILE_ACCESS);
2582 H5Pset_fapl_mpio(mappings_plist_id, PROTEUS_COMM_WORLD, MPI_INFO_NULL);
2583 file_id = H5Fopen(H5FILE_NAME, H5F_ACC_RDWR, mappings_plist_id);
2584 H5Pclose(mappings_plist_id);
2585 hid_t eb_new2old_dataset_id = H5Dcreate(file_id,
"elementBoundaryNumbering_new2old", H5T_NATIVE_INT,
2586 eb_new2old_filespace_id,
2587 H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
2588 hsize_t eb_count[] = {
static_cast<hsize_t
>(nElementBoundaries_subdomain_new[rank])};
2589 hsize_t eb_offset[] = {
static_cast<hsize_t
>(elementBoundaryOffsets_new[rank])};
2590 hid_t eb_new2old_memspace_id = H5Screate_simple(ARRAY_RANK, eb_count, NULL);
2591 H5Sselect_hyperslab(eb_new2old_filespace_id, H5S_SELECT_SET, eb_offset, NULL, eb_count, NULL);
2592 hid_t eb_new2old_plist_id = H5Pcreate(H5P_DATASET_XFER);
2593 H5Pset_dxpl_mpio(eb_new2old_plist_id, H5FD_MPIO_COLLECTIVE);
2594 status = H5Dwrite(eb_new2old_dataset_id, H5T_NATIVE_INT, eb_new2old_memspace_id, eb_new2old_filespace_id,
2595 eb_new2old_plist_id, &elementBoundaryNumbering_new2old[0]);
2596 H5Pclose(eb_new2old_plist_id);
2597 H5Sclose(eb_new2old_memspace_id);
2598 H5Sclose(eb_new2old_filespace_id);
2599 H5Dclose(eb_new2old_dataset_id);
2601 valarray<hsize_t> old_elementBoundary_indices(nElementBoundaries_subdomain_new[rank]);
2602 valarray<int> new_elementBoundary_indices(nElementBoundaries_subdomain_new[rank]);
2603 for (
int i=0;i<nElementBoundaries_subdomain_new[rank];i++)
2605 old_elementBoundary_indices[i] = elementBoundaryNumbering_new2old[i];
2606 new_elementBoundary_indices[i] = elementBoundaryOffsets_new[rank]+i;
2608 hid_t eb_old2new_filespace_id = H5Screate_simple(ARRAY_RANK, eb_dims, NULL);
2609 hid_t eb_old2new_dataset_id = H5Dcreate(file_id,
"elementBoundaryNumbering_old2new", H5T_NATIVE_INT,
2610 eb_old2new_filespace_id,
2611 H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
2612 hid_t eb_old2new_memspace_id = H5Screate_simple(ARRAY_RANK, eb_count, NULL);
2613 status = H5Sselect_elements(eb_old2new_filespace_id, H5S_SELECT_SET,
2614 nElementBoundaries_subdomain_new[rank], &old_elementBoundary_indices[0]);
2615 hid_t eb_old2new_plist_id = H5Pcreate(H5P_DATASET_XFER);
2616 H5Pset_dxpl_mpio(eb_old2new_plist_id, H5FD_MPIO_COLLECTIVE);
2617 status = H5Dwrite(eb_old2new_dataset_id, H5T_NATIVE_INT,
2618 eb_old2new_memspace_id, eb_old2new_filespace_id,
2619 eb_old2new_plist_id, &new_elementBoundary_indices[0]);
2620 H5Pclose(eb_new2old_plist_id);
2621 H5Sclose(eb_old2new_memspace_id);
2622 H5Dclose(eb_old2new_dataset_id);
2623 H5Sclose(eb_old2new_filespace_id);
2627 valarray<hsize_t> old_elementBoundary_indices_subdomain(elementBoundaries_subdomain.size());
2628 valarray<int> new_elementBoundary_indices_subdomain(elementBoundaries_subdomain.size());
2629 int ebN_subdomain = 0;
2630 for (
auto it = elementBoundaries_subdomain.begin();
2631 it != elementBoundaries_subdomain.end();
2632 it++, ebN_subdomain++)
2634 old_elementBoundary_indices_subdomain[ebN_subdomain] =
static_cast<hsize_t
>(*it);
2636 file_id = H5Fopen(H5FILE_NAME, H5F_ACC_RDONLY, H5P_DEFAULT);
2637 eb_old2new_filespace_id = H5Screate_simple(ARRAY_RANK, eb_dims, NULL);
2638 eb_old2new_dataset_id = H5Dopen1(file_id,
"/elementBoundaryNumbering_old2new");
2639 status = H5Sselect_elements(eb_old2new_filespace_id, H5S_SELECT_SET,
2640 elementBoundaries_subdomain.size(),
2641 &old_elementBoundary_indices_subdomain[0]);
2642 hsize_t eb_subdomain_count[]={
static_cast<hsize_t
>(elementBoundaries_subdomain.size())};
2643 hid_t eb_old2new_subdomain_memspace_id = H5Screate_simple(ARRAY_RANK, eb_subdomain_count, NULL);
2644 hid_t eb_old2new_subdomain_plist_id = H5Pcreate(H5P_DATASET_XFER);
2645 H5Pset_dxpl_mpio(eb_old2new_subdomain_plist_id, H5FD_MPIO_COLLECTIVE);
2646 status = H5Dread(eb_old2new_dataset_id, H5T_NATIVE_INT,
2647 eb_old2new_subdomain_memspace_id, eb_old2new_filespace_id,
2648 H5P_DEFAULT, &new_elementBoundary_indices_subdomain[0]);
2649 H5Pclose(eb_old2new_subdomain_plist_id);
2650 H5Sclose(eb_old2new_subdomain_memspace_id);
2651 H5Sclose(eb_old2new_filespace_id);
2652 H5Dclose(eb_old2new_dataset_id);
2654 map<int,int> elementBoundaryNumbering_old2new_subdomain_map;
2655 map<int,int> elementBoundaryNumbering_new2old_subdomain_map;
2656 for (
int i=0;i<elementBoundaries_subdomain.size();i++)
2658 elementBoundaryNumbering_old2new_subdomain_map[old_elementBoundary_indices_subdomain[i]] = new_elementBoundary_indices_subdomain[i];
2659 elementBoundaryNumbering_new2old_subdomain_map[new_elementBoundary_indices_subdomain[i]] = old_elementBoundary_indices_subdomain[i];
2661 IS elementBoundaryNumberingIS_subdomain_new2old;
2662 ISCreateGeneral(PROTEUS_COMM_WORLD,elementBoundaries_subdomain_owned.size(),&elementBoundaryNumbering_new2old[0],PETSC_COPY_VALUES,&elementBoundaryNumberingIS_subdomain_new2old);
2665 IS elementBoundaryNumberingIS_global_new2old;
2666 const PetscInt *elementBoundaryNumbering_global_new2old;
2667 valarray<int> elementBoundaryNumbering_global_old2new(nElementBoundaries_global);
2668 ISAllGather(elementBoundaryNumberingIS_subdomain_new2old,&elementBoundaryNumberingIS_global_new2old);
2670 ierr =
enforceMemoryLimit(PROTEUS_COMM_WORLD, rank, max_rss_gb,
"Allocating elementBoudnary old2new/new2old");CHKERRABORT(PROTEUS_COMM_WORLD, ierr);
2671 ISGetIndices(elementBoundaryNumberingIS_global_new2old,&elementBoundaryNumbering_global_new2old);
2672 for (
int i=0;i<nElementBoundaries_global;i++)
2674 elementBoundaryNumbering_global_old2new[elementBoundaryNumbering_global_new2old[i]]=i;
2678 valarray<int> elementBoundaryNumbering_old2new_read(nElementBoundaries_global);
2679 file_id = H5Fopen(H5FILE_NAME, H5F_ACC_RDONLY, H5P_DEFAULT);
2680 eb_old2new_dataset_id = H5Dopen1(file_id,
"/elementBoundaryNumbering_old2new");
2681 status = H5Dread(eb_old2new_dataset_id, H5T_NATIVE_INT, H5S_ALL, H5S_ALL, H5P_DEFAULT,
2682 &elementBoundaryNumbering_old2new_read[0]);
2683 status = H5Dclose(eb_old2new_dataset_id);
2684 for (
int i=0;i<nElementBoundaries_global;i++)
2685 assert(elementBoundaryNumbering_global_old2new[i] == elementBoundaryNumbering_old2new_read[i]);
2686 std::cout<<
"==================out of core old2new elementBoundaries is correct!===================="<<std::endl;
2687 hid_t eb_new2old_dataset_id = H5Dopen1(file_id,
"/elementBoundaryNumbering_new2old");
2688 valarray<int> elementBoundaryNumbering_new2old_read(nElementBoundaries_global);
2689 status = H5Dread(eb_new2old_dataset_id, H5T_NATIVE_INT, H5S_ALL, H5S_ALL, H5P_DEFAULT,
2690 &elementBoundaryNumbering_new2old_read[0]);
2691 status = H5Dclose(eb_new2old_dataset_id);
2692 status = H5Fclose(file_id);
2693 for (
int i=0;i<nElementBoundaries_global;i++)
2694 assert(elementBoundaryNumbering_global_new2old[i] == elementBoundaryNumbering_new2old_read[i]);
2695 std::cout<<
"==================out of core new2old elementBoundaries is correct!===================="<<std::endl;
2696 for (
int i=0;i<elementBoundaries_subdomain.size();i++)
2698 int ebN = old_elementBoundary_indices_subdomain[i];
2699 assert(elementBoundaryNumbering_global_old2new[ebN] == elementBoundaryNumbering_old2new_subdomain_map[ebN]);
2701 std::cout<<
"==================out of core elementBoundaries old2new subdomain_map is correct!===================="<<std::endl;
2702 for (
int i=0;i<elementBoundaries_subdomain.size();i++)
2704 int ebN = new_elementBoundary_indices_subdomain[i];
2705 assert(elementBoundaryNumbering_global_new2old[ebN] == elementBoundaryNumbering_new2old_subdomain_map[ebN]);
2707 std::cout<<
"==================out of core elementBoundaries new2old subdomain_map is correct!===================="<<std::endl;
2708 ISRestoreIndices(elementBoundaryNumberingIS_global_new2old,&elementBoundaryNumbering_global_new2old);
2709 ierr=ISDestroy(&elementBoundaryNumberingIS_subdomain_new2old);CHKERRABORT(PROTEUS_COMM_WORLD, ierr);
2710 ierr=ISDestroy(&elementBoundaryNumberingIS_global_new2old);CHKERRABORT(PROTEUS_COMM_WORLD, ierr);
2712 PetscLogEventEnd(build_subdomains_faces_event,0,0,0,0);
2713 ierr =
enforceMemoryLimit(PROTEUS_COMM_WORLD, rank, max_rss_gb,
"Done allocating elementBoudnary old2new/new2old");CHKERRABORT(PROTEUS_COMM_WORLD, ierr);
2714 int build_subdomains_edges_event;
2715 PetscLogEventRegister(
"Subd edges",0,&build_subdomains_edges_event);
2716 PetscLogEventBegin(build_subdomains_edges_event,0,0,0,0);
2721 logEvent(
"Generating global edge numbering and ownership corresponding to new subdomain ownership",5);
2722 std::ifstream edgeFile(edgeFileName.c_str());
2724 if (!edgeFile.good())
2726 std::cerr<<
"cannot open Tetgen edge file"
2727 <<edgeFileName<<std::endl;
2732 bool hasEdgeMarkers =
false;
2734 int ihasEdgeMarkers(0);
2736 assert(nEdges_global > 0);
2737 if (ihasEdgeMarkers > 0)
2739 hasEdgeMarkers =
true;
2742 set<int> edges_subdomain_owned;
2744 map<int,int> edgeMaterialTypesMap;
2745 map<int,valarray<int> > elementEdgesMap;
2746 map<int,pair<int,int> > edgeNodesMap;
2747 set<int> supportedEdges;
2748 for (
int ied = 0; ied < nEdges_global; ied++)
2750 int ned,nn0,nn1;
int edId(0);
2752 if (ihasEdgeMarkers > 0)
2757 assert(0 <= ned && ned < nEdges_global && edgeFile.good());
2759 if (nodes_old2new_subdomain_map.find(nn0) != nodes_old2new_subdomain_map.end())
2760 nn0_new = nodes_old2new_subdomain_map[nn0];
2761 if (nn0_new >= nodeOffsets_new[rank] && nn0_new < nodeOffsets_new[rank+1])
2763 nodeEdgesStar.at(nn0_new-nodeOffsets_new[rank]).insert(ned);
2764 supportedEdges.insert(ned);
2767 if(nodes_old2new_subdomain_map.find(nn1) != nodes_old2new_subdomain_map.end())
2768 nn1_new = nodes_old2new_subdomain_map[nn1];
2769 if (nn1_new >= nodeOffsets_new[rank] && nn1_new < nodeOffsets_new[rank+1])
2771 nodeEdgesStar.at(nn1_new-nodeOffsets_new[rank]).insert(ned);
2772 supportedEdges.insert(ned);
2774 int nodes[2] = {nn0_new,nn1_new};
2777 if (edgeElementsMap.find(nodeTuple) != edgeElementsMap.end())
2779 if (nodeTuple.
nodes[0] >= nodeOffsets_new[rank] && nodeTuple.
nodes[0] < nodeOffsets_new[rank+1])
2780 edges_subdomain_owned.insert(ned);
2782 edgeNodesMap[ned].first = nodeTuple.
nodes[0];
2783 edgeNodesMap[ned].second = nodeTuple.
nodes[1];
2784 if (ihasEdgeMarkers > 0)
2785 edgeMaterialTypesMap[ned]=edId;
2786 for (
auto elementp=edgeElementsMap[nodeTuple].begin();
2787 elementp != edgeElementsMap[nodeTuple].end();
2790 int eN = elementNumbering_old2new_subdomain_map[elementp->first];
2791 if (elementEdgesMap.find(eN) != elementEdgesMap.end())
2793 elementEdgesMap[eN][elementp->second] = ned;
2797 std::valarray<int>
init(-1,6);
2798 elementEdgesMap[eN] =
init;
2799 elementEdgesMap[eN][elementp->second] = ned;
2805 int nEdges_owned_subdomain=edges_subdomain_owned.size(),
2807 MPI_Allreduce(&nEdges_owned_subdomain,&nEdges_owned_new,1,MPI_INT,MPI_SUM,PROTEUS_COMM_WORLD);
2808 assert(nEdges_owned_new == nEdges_global);
2809 ierr =
enforceMemoryLimit(PROTEUS_COMM_WORLD, rank, max_rss_gb,
"Done reading edges");CHKERRABORT(PROTEUS_COMM_WORLD, ierr);
2815 for (
auto edgesp=elementEdgesMap.begin();
2816 edgesp!=elementEdgesMap.end();
2820 for (
int iv=0;iv<4;iv++)
2824 if (elementNumbering_new2old_subdomain_map.find(edgesp->first) != elementNumbering_new2old_subdomain_map.end())
2825 nN_global = elementNodesArrayMap[elementNumbering_new2old_subdomain_map[edgesp->first]][iv];
2826 if (nN_global >= nodeOffsets_new[rank] && nN_global < nodeOffsets_new[rank+1])
2829 for(
int ed=0;ed<6;ed++)
2831 nodeEdgesStar.at(nN_global-nodeOffsets_new[rank]).insert(edgesp->second[ed]);
2837 valarray<int> nodeEdgeOffsets(nNodes_subdomain_new[rank]+1);
2838 nodeEdgeOffsets[0] = 0;
2839 for (
int nN = 0; nN < nNodes_subdomain_new[rank]; nN++)
2840 nodeEdgeOffsets[nN+1] = nodeEdgeOffsets[nN]+nodeEdgesStar.at(nN).size();
2841 valarray<int> nodeEdgesArray(nodeEdgeOffsets[nNodes_subdomain_new[rank]]);
2842 for (
int nN=0,offset=0; nN < nNodes_subdomain_new[rank]; nN++)
2844 for (
auto edN_star = nodeEdgesStar.at(nN).begin();
2845 edN_star != nodeEdgesStar.at(nN).end();
2846 edN_star++,offset++)
2848 nodeEdgesArray[offset] = *edN_star;
2853 valarray<int> nEdges_subdomain_new(size),
2854 edgeOffsets_new(size+1);
2855 for (
int sdN=0;sdN<size;sdN++)
2857 nEdges_subdomain_new[sdN] = edges_subdomain_owned.size();
2859 nEdges_subdomain_new[sdN] = 0;
2860 valarray<int> nEdges_subdomain_new_send=nEdges_subdomain_new;
2861 MPI_Allreduce(&nEdges_subdomain_new_send[0],&nEdges_subdomain_new[0],size,MPI_INT,MPI_SUM,PROTEUS_COMM_WORLD);
2865 edgeOffsets_new[0] = 0;
2866 for (
int sdN=0;sdN<size;sdN++)
2867 edgeOffsets_new[sdN+1] = edgeOffsets_new[sdN]+nEdges_subdomain_new[sdN];
2873 valarray<int> edgeNumbering_new2old(edges_subdomain_owned.size());
2874 map<int,int> edgeNumbering_old2new_map;
2875 set<int>::iterator edN_ownedp=edges_subdomain_owned.begin();
2876 for (
int edN=0;edN<int(edges_subdomain_owned.size());edN++,edN_ownedp++)
2878 edgeNumbering_new2old[edN] = *edN_ownedp;
2879 edgeNumbering_old2new_map[*edN_ownedp] = edN+edgeOffsets_new[rank];
2884 logEvent(
"Writing/reading edge numberings to hdf5",5);
2885 hsize_t ed_dims[]={
static_cast<hsize_t
>(nEdges_global)};
2886 hid_t ed_new2old_filespace_id = H5Screate_simple(ARRAY_RANK, ed_dims, NULL);
2887 mappings_plist_id = H5Pcreate(H5P_FILE_ACCESS);
2888 H5Pset_fapl_mpio(mappings_plist_id, PROTEUS_COMM_WORLD, MPI_INFO_NULL);
2889 file_id = H5Fopen(H5FILE_NAME, H5F_ACC_RDWR, mappings_plist_id);
2890 H5Pclose(mappings_plist_id);
2891 hid_t ed_new2old_dataset_id = H5Dcreate(file_id,
"edgeNumbering_new2old", H5T_NATIVE_INT,
2892 ed_new2old_filespace_id,
2893 H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
2894 hsize_t ed_count[]={
static_cast<hsize_t
>(nEdges_subdomain_new[rank])};
2895 hsize_t ed_offset[]={
static_cast<hsize_t
>(edgeOffsets_new[rank])};
2896 hid_t ed_new2old_memspace_id = H5Screate_simple(ARRAY_RANK, ed_count, NULL);
2897 H5Sselect_hyperslab(ed_new2old_filespace_id, H5S_SELECT_SET, ed_offset, NULL, ed_count, NULL);
2898 hid_t ed_new2old_plist_id = H5Pcreate(H5P_DATASET_XFER);
2899 H5Pset_dxpl_mpio(ed_new2old_plist_id, H5FD_MPIO_COLLECTIVE);
2900 status = H5Dwrite(ed_new2old_dataset_id, H5T_NATIVE_INT,
2901 ed_new2old_memspace_id, ed_new2old_filespace_id,
2902 ed_new2old_plist_id, &edgeNumbering_new2old[0]);
2903 H5Pclose(ed_new2old_plist_id);
2904 H5Sclose(ed_new2old_memspace_id);
2905 H5Sclose(ed_new2old_filespace_id);
2906 H5Dclose(ed_new2old_dataset_id);
2908 valarray<hsize_t> old_edge_indices(
static_cast<hsize_t
>(nEdges_subdomain_new[rank]));
2909 valarray<int> new_edge_indices(nEdges_subdomain_new[rank]);
2910 for (
int i=0;i<nEdges_subdomain_new[rank];i++)
2912 old_edge_indices[i] =
static_cast<hsize_t
>(edgeNumbering_new2old[i]);
2913 new_edge_indices[i] = edgeOffsets_new[rank]+i;
2915 hid_t ed_old2new_filespace_id = H5Screate_simple(ARRAY_RANK, ed_dims, NULL);
2916 hid_t ed_old2new_dataspace_id = H5Dcreate(file_id,
"edgeNumbering_old2new", H5T_NATIVE_INT,
2917 ed_old2new_filespace_id,
2918 H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
2919 hid_t ed_old2new_memspace_id = H5Screate_simple(ARRAY_RANK, ed_count, NULL);
2920 ed_old2new_filespace_id = H5Dget_space(ed_old2new_dataspace_id);
2921 status = H5Sselect_elements(ed_old2new_filespace_id, H5S_SELECT_SET,
2922 nEdges_subdomain_new[rank], &old_edge_indices[0]);
2923 hid_t ed_old2new_plist_id = H5Pcreate(H5P_DATASET_XFER);
2924 H5Pset_dxpl_mpio(ed_old2new_plist_id, H5FD_MPIO_COLLECTIVE);
2925 status = H5Dwrite(ed_old2new_dataspace_id, H5T_NATIVE_INT,
2926 ed_old2new_memspace_id, ed_old2new_filespace_id,
2927 ed_old2new_plist_id, &new_edge_indices[0]);
2928 H5Pclose(ed_old2new_plist_id);
2929 H5Sclose(ed_old2new_memspace_id);
2930 H5Dclose(ed_old2new_dataspace_id);
2931 H5Sclose(ed_old2new_filespace_id);
2935 valarray<hsize_t> old_edge_indices_subdomain(edgeNodesMap.size());
2936 valarray<int> new_edge_indices_subdomain(edgeNodesMap.size());
2937 int edN_old_subdomain = 0;
2938 for (
auto it = edgeNodesMap.begin(); it != edgeNodesMap.end(); it++, edN_old_subdomain++)
2940 old_edge_indices_subdomain[edN_old_subdomain] =
static_cast<hsize_t
>(it->first);
2942 file_id = H5Fopen(H5FILE_NAME, H5F_ACC_RDONLY, H5P_DEFAULT);
2943 ed_old2new_dataspace_id = H5Dopen1(file_id,
"/edgeNumbering_old2new");
2944 ed_old2new_filespace_id = H5Screate_simple(ARRAY_RANK, ed_dims, NULL);
2945 status = H5Sselect_elements(ed_old2new_filespace_id, H5S_SELECT_SET,
2946 edgeNodesMap.size(), &old_edge_indices_subdomain[0]);
2947 hsize_t ed_subdomain_count[]={
static_cast<hsize_t
>(edgeNodesMap.size())};
2948 hid_t ed_old2new_subdomain_memspace_id = H5Screate_simple(ARRAY_RANK, ed_subdomain_count, NULL);
2949 hid_t ed_old2new_subdomain_plist_id = H5Pcreate(H5P_DATASET_XFER);
2950 H5Pset_dxpl_mpio(ed_old2new_subdomain_plist_id, H5FD_MPIO_COLLECTIVE);
2951 status = H5Dread(ed_old2new_dataspace_id, H5T_NATIVE_INT, ed_old2new_subdomain_memspace_id,
2952 ed_old2new_filespace_id,
2953 H5P_DEFAULT, &new_edge_indices_subdomain[0]);
2954 H5Pclose(ed_old2new_subdomain_plist_id);
2955 H5Sclose(ed_old2new_subdomain_memspace_id);
2956 H5Sclose(ed_old2new_filespace_id);
2957 H5Dclose(ed_old2new_dataspace_id);
2959 map<int,int> edgeNumbering_old2new_subdomain_map;
2960 map<int,int> edgeNumbering_new2old_subdomain_map;
2961 for (
int i=0;i<edgeNodesMap.size();i++)
2963 edgeNumbering_old2new_subdomain_map[old_edge_indices_subdomain[i]] = new_edge_indices_subdomain[i];
2964 edgeNumbering_new2old_subdomain_map[new_edge_indices_subdomain[i]] = old_edge_indices_subdomain[i];
2968 IS edgeNumberingIS_subdomain_new2old;
2969 ISCreateGeneral(PROTEUS_COMM_WORLD,edges_subdomain_owned.size(),&edgeNumbering_new2old[0],PETSC_COPY_VALUES,&edgeNumberingIS_subdomain_new2old);
2970 IS edgeNumberingIS_global_new2old;
2971 ISAllGather(edgeNumberingIS_subdomain_new2old,&edgeNumberingIS_global_new2old);
2972 const PetscInt *edgeNumbering_global_new2old;
2973 ISGetIndices(edgeNumberingIS_global_new2old,&edgeNumbering_global_new2old);
2974 valarray<int> edgeNumbering_global_old2new(newMesh.
nEdges_global);
2975 ierr =
enforceMemoryLimit(PROTEUS_COMM_WORLD, rank, max_rss_gb,
"Setting edgeNumering old2new/new2old");CHKERRABORT(PROTEUS_COMM_WORLD, ierr);
2978 edgeNumbering_global_old2new[edgeNumbering_global_new2old[edN]] = edN;
2982 valarray<int> edgeNumbering_old2new_read(nEdges_global);
2983 file_id = H5Fopen(H5FILE_NAME, H5F_ACC_RDONLY, H5P_DEFAULT);
2984 hid_t ed_old2new_dataset_id = H5Dopen1(file_id,
"/edgeNumbering_old2new");
2985 status = H5Dread(ed_old2new_dataset_id, H5T_NATIVE_INT, H5S_ALL, H5S_ALL, H5P_DEFAULT,
2986 &edgeNumbering_old2new_read[0]);
2987 status = H5Dclose(ed_old2new_dataset_id);
2988 for (
int i=0;i<nEdges_global;i++)
2989 assert(edgeNumbering_global_old2new[i] == edgeNumbering_old2new_read[i]);
2990 std::cout<<
"==================out of core old2new edges is correct!===================="<<std::endl;
2991 hid_t ed_new2old_dataset_id = H5Dopen1(file_id,
"/edgeNumbering_new2old");
2992 valarray<int> edgeNumbering_new2old_read(nEdges_global);
2993 status = H5Dread(ed_new2old_dataset_id, H5T_NATIVE_INT, H5S_ALL, H5S_ALL, H5P_DEFAULT,
2994 &edgeNumbering_new2old_read[0]);
2995 status = H5Dclose(ed_new2old_dataset_id);
2996 status = H5Fclose(file_id);
2997 for (
int i=0;i<nEdges_global;i++)
2998 assert(edgeNumbering_global_new2old[i] == edgeNumbering_new2old_read[i]);
2999 std::cout<<
"==================out of core new2old edges is correct!===================="<<std::endl;
3001 for (
auto it = edgeNodesMap.begin(); it != edgeNodesMap.end(); ++it)
3003 assert(edgeNumbering_global_old2new[it->first] == edgeNumbering_old2new_subdomain_map[it->first]);
3005 std::cout<<
"==================out of core edgeNumbering old2new subdomain is correct!===================="<<std::endl;
3006 for (
int i = 0; i < edgeNodesMap.size(); i++)
3008 int edN = new_edge_indices_subdomain[i];
3009 assert(edgeNumbering_global_new2old[edN] == edgeNumbering_new2old_subdomain_map[edN]);
3011 std::cout<<
"==================out of core edgeNumbering new2old subdomain is correct!===================="<<std::endl;
3012 ISRestoreIndices(edgeNumberingIS_global_new2old,&edgeNumbering_global_new2old);
3013 ierr=ISDestroy(&edgeNumberingIS_subdomain_new2old);CHKERRABORT(PROTEUS_COMM_WORLD, ierr);
3014 ierr=ISDestroy(&edgeNumberingIS_global_new2old);CHKERRABORT(PROTEUS_COMM_WORLD, ierr);
3016 ierr =
enforceMemoryLimit(PROTEUS_COMM_WORLD, rank, max_rss_gb,
"Done allocating edgeNumering old2new/new2old");CHKERRABORT(PROTEUS_COMM_WORLD, ierr);
3020 logEvent(
"Creating ghost information",5);
3021 set<int> elements_overlap,nodes_overlap,elementBoundaries_overlap,edges_overlap;
3022 for (
int nN = 0; nN < nNodes_subdomain_new[rank]; nN++)
3025 for (
int offset = nodeStarOffsetsNew[nN];offset<nodeStarOffsetsNew[nN+1];offset++)
3027 int nN_point_global = nodeStarArrayNew[offset];
3028 bool offproc = nN_point_global < nodeOffsets_new[rank] || nN_point_global >= nodeOffsets_new[rank+1];
3030 nodes_overlap.insert(nN_point_global);
3033 for (
int eN_star_offset = nodeElementOffsets[nN];
3034 eN_star_offset < nodeElementOffsets[nN+1]; eN_star_offset++)
3036 int eN_star_old = nodeElementsArray[eN_star_offset];
3037 int eN_star_new = elementNumbering_old2new_subdomain_map[eN_star_old];
3038 bool offproc = eN_star_new >= elementOffsets_new[rank+1] || eN_star_new < elementOffsets_new[rank];
3040 elements_overlap.insert(eN_star_new);
3043 for (
int ebN_star_offset = nodeElementBoundaryOffsets[nN];
3044 ebN_star_offset < nodeElementBoundaryOffsets[nN+1]; ebN_star_offset++)
3046 int ebN_star_old = nodeElementBoundariesArray[ebN_star_offset];
3047 int ebN_star_new = elementBoundaryNumbering_old2new_subdomain_map[ebN_star_old];
3048 bool offproc = ebN_star_new >= elementBoundaryOffsets_new[rank+1] || ebN_star_new < elementBoundaryOffsets_new[rank];
3050 elementBoundaries_overlap.insert(ebN_star_new);
3053 for (
int edN_star_offset = nodeEdgeOffsets[nN];
3054 edN_star_offset < nodeEdgeOffsets[nN+1]; edN_star_offset++)
3056 int edN_star_old = nodeEdgesArray[edN_star_offset];
3057 int edN_star_new = edgeNumbering_old2new_subdomain_map[edN_star_old];
3058 bool offproc = edN_star_new >= edgeOffsets_new[rank+1] || edN_star_new < edgeOffsets_new[rank];
3060 edges_overlap.insert(edN_star_new);
3064 assert(edges_overlap.size() + nEdges_subdomain_new[rank] == edgeNodesMap.size());
3068 int nN_subdomain = nNodes_subdomain_new[rank];
3069 map<int,int> nodes_overlap_global2subdomainMap;
3070 for (
auto nN_globalp=nodes_overlap.begin();nN_globalp != nodes_overlap.end(); nN_globalp++,nN_subdomain++)
3071 nodes_overlap_global2subdomainMap[*nN_globalp] = nN_subdomain;
3076 assert(nNodes_overlap<2);
3077 logEvent(
"Skipping additiona of addtional overlap--only minimal overlap supported by this partioning function for now",5);
3079 PetscLogEventEnd(build_subdomains_edges_event,0,0,0,0);
3080 int build_subdomains_renumber_event;
3081 PetscLogEventRegister(
"Subd's renumber",0,&build_subdomains_renumber_event);
3082 PetscLogEventBegin(build_subdomains_renumber_event,0,0,0,0);
3086 logEvent(
"USER WARNING: You are partioning in parallel directly from Tetgen files (parun -F ...). ",5);
3087 logEvent(
"USER WARNING: You must use the 'f' and 'ee' flags in the triangleOptions input to avoid errors.",5);
3088 logEvent(
"USER WARNING: You may have done so, but the tetgen options are not stored in the Tetgen files, so we can't know.",5);
3089 logEvent(
"Creating final subdomain mesh data structures",5);
3105 map<int,int> nodeNumbering_global2subdomainMap;
3106 map<int,int> elementBoundaryNumbering_global2subdomainMap;
3107 map<int,int> edgeNumbering_global2subdomainMap;
3113 for (
int iv = 0; iv < nNodes_global; iv++)
3115 int nv;
double x,y,
z;
int nodeId(0);
3117 if (hasVertexMarkers > 0)
3118 vertexFile >> nodeId;
3120 assert(0 <= nv && nv < nNodes_global && vertexFile.good());
3121 int nN_global_new = -1;
3122 if (nodes_old2new_subdomain_map.find(nv) != nodes_old2new_subdomain_map.end())
3123 nN_global_new = nodes_old2new_subdomain_map[nv];
3125 if (nN_global_new >= nodeOffsets_new[rank] && nN_global_new < nodeOffsets_new[rank+1])
3127 int nv_subdomain_new = nN_global_new - nodeOffsets_new[rank];
3128 nodeNumbering_subdomain2global[nv_subdomain_new] = nN_global_new;
3129 nodeNumbering_global2subdomainMap[nN_global_new] = nv_subdomain_new;
3133 if (hasVertexMarkers > 0)
3137 if (nodes_overlap.count(nN_global_new) == 1)
3139 int nv_subdomain_new = nodes_overlap_global2subdomainMap[nN_global_new];
3140 nodeNumbering_subdomain2global[nv_subdomain_new] = nN_global_new;
3141 nodeNumbering_global2subdomainMap[nN_global_new] = nv_subdomain_new;
3145 if (hasVertexMarkers > 0)
3152 ierr =
enforceMemoryLimit(PROTEUS_COMM_WORLD, rank, max_rss_gb,
"Done reading vertices");CHKERRABORT(PROTEUS_COMM_WORLD, ierr);
3161 for (
int eN = 0; eN < nElements_subdomain_new[rank]; eN++)
3163 int eN_global_new = elementOffsets_new[rank] + eN;
3164 int eN_global_old = elementNumbering_subdomain_new2old[eN];
3165 elementNumbering_subdomain2global[eN] = eN_global_new;
3169 int nN_global_new = elementNodesArrayMap[eN_global_old][nN];
3170 int nN_subdomain = nodeNumbering_global2subdomainMap[nN_global_new];
3177 set<int>::iterator eN_p = elements_overlap.begin();
3178 for (
int eN = nElements_subdomain_new[rank]; eN < nElements_subdomain_new[rank] + int(elements_overlap.size()); eN++,eN_p++)
3180 int eN_global_new = *eN_p;
3181 int eN_global_old = elementNumbering_new2old_subdomain_map[eN_global_new];
3182 elementNumbering_subdomain2global[eN] = eN_global_new;
3186 int nN_global_new = elementNodesArrayMap[eN_global_old][nN];
3187 int nN_subdomain = nodeNumbering_global2subdomainMap[nN_global_new];
3196 for (
int ebN=0; ebN < nElementBoundaries_subdomain_new[rank]; ebN++)
3198 int ebN_global = ebN + elementBoundaryOffsets_new[rank];
3199 elementBoundaryNumbering_subdomain2global[ebN]=ebN_global;
3200 elementBoundaryNumbering_global2subdomainMap[ebN_global] = ebN;
3205 set<int>::iterator ebN_p = elementBoundaries_overlap.begin();
3206 for(
int ebN=nElementBoundaries_subdomain_new[rank];ebN < nElementBoundaries_subdomain_new[rank] + int(elementBoundaries_overlap.size()); ebN++,ebN_p++)
3208 int ebN_global = *ebN_p;
3209 elementBoundaryNumbering_subdomain2global[ebN] = ebN_global;
3210 elementBoundaryNumbering_global2subdomainMap[ebN_global] = ebN;
3219 for (
int eN=0;eN<nElements_subdomain_new[rank];eN++)
3221 int eN_global = eN+elementOffsets_new[rank];
3225 elementBoundaryNumbering_global2subdomainMap[elementBoundaryNumbering_old2new_subdomain_map[elementBoundariesMap[eN_global][ebN]]];
3231 eN_p = elements_overlap.begin();
3232 for (
int eN = nElements_subdomain_new[rank]; eN < nElements_subdomain_new[rank] + int(elements_overlap.size()); eN++,eN_p++)
3234 int eN_global_new = *eN_p;
3238 elementBoundaryNumbering_global2subdomainMap[elementBoundaryNumbering_old2new_subdomain_map[elementBoundariesMap[eN_global_new][ebN]]];
3246 for (
int edN=0; edN < nEdges_subdomain_new[rank]; edN++)
3248 int edN_global = edN + edgeOffsets_new[rank];
3249 edgeNumbering_subdomain2global[edN]=edN_global;
3250 edgeNumbering_global2subdomainMap[edN_global] = edN;
3255 set<int>::iterator edN_p = edges_overlap.begin();
3256 for(
int edN=nEdges_subdomain_new[rank];edN < nEdges_subdomain_new[rank] + int(edges_overlap.size()); edN++,edN_p++)
3258 int edN_global = *edN_p;
3259 edgeNumbering_subdomain2global[edN] = edN_global;
3260 edgeNumbering_global2subdomainMap[edN_global] = edN;
3268 for (map<
int,pair<int,int> >::iterator edgep=edgeNodesMap.begin();
3269 edgep!=edgeNodesMap.end();
3272 int edN_global_old = edgep->first;
3273 int edN_global_new = edgeNumbering_old2new_subdomain_map[edN_global_old];
3274 assert(edgeNumbering_global2subdomainMap.find(edN_global_new) != edgeNumbering_global2subdomainMap.end());
3275 int edN_subdomain = edgeNumbering_global2subdomainMap[edN_global_new];
3288 if (hasElementBoundaryMarkers)
3291 for (
auto ebmp = elementBoundaryMaterialTypesMap.begin(); ebmp != elementBoundaryMaterialTypesMap.end();ebmp++)
3293 int ebN_global_new = elementBoundaryNumbering_old2new_subdomain_map[ebmp->first];
3294 assert(elementBoundaryNumbering_global2subdomainMap.find(ebN_global_new) != elementBoundaryNumbering_global2subdomainMap.end());
3295 int ebN_subdomain = elementBoundaryNumbering_global2subdomainMap[ebN_global_new];
3298 if (!hasVertexMarkers)
3309 ierr =
enforceMemoryLimit(PROTEUS_COMM_WORLD, rank, max_rss_gb,
"Done with material types");CHKERRABORT(PROTEUS_COMM_WORLD, ierr);
3310 PetscLogEventEnd(build_subdomains_renumber_event,0,0,0,0);
3311 int build_subdomains_cleanup_event;
3312 PetscLogEventRegister(
"Cleanup",0,&build_subdomains_cleanup_event);
3313 PetscLogEventBegin(build_subdomains_cleanup_event,0,0,0,0);
3315 logEvent(
"Cleaning up after partitioning",5);
3328 for (
int sdN = 0; sdN < size+1; sdN++)
3358 PetscLogEventEnd(build_subdomains_cleanup_event,0,0,0,0);
3360 PetscLogView(PETSC_VIEWER_STDOUT_WORLD);
3361 ierr =
enforceMemoryLimit(PROTEUS_COMM_WORLD, rank, max_rss_gb,
"Done with partitioning!");CHKERRABORT(PROTEUS_COMM_WORLD, ierr);
3366 using namespace std;
3367 PetscErrorCode ierr;
3368 PetscMPIInt size,rank;
3370 ierr = MPI_Comm_size(PROTEUS_COMM_WORLD,&size);CHKERRABORT(PROTEUS_COMM_WORLD, ierr);
3371 ierr = MPI_Comm_rank(PROTEUS_COMM_WORLD,&rank);CHKERRABORT(PROTEUS_COMM_WORLD, ierr);
3372 PetscLogStage partitioning_stage;
3373 PetscLogStageRegister(
"Mesh Partition",&partitioning_stage);
3374 PetscLogStagePush(partitioning_stage);
3401 bool failed =
false;
3402 const int simplexDim = 3;
3403 const int vertexDim = 3;
3405 std::string vertexFileName = std::string(filebase) +
".node" ;
3406 std::string elementFileName = std::string(filebase) +
".ele" ;
3407 std::string elementBoundaryFileName = std::string(filebase) +
".edge" ;
3408 std::string edgeFileName = std::string(filebase) +
".edge" ;
3419 int read_elements_event;
3420 PetscLogEventRegister(
"Read eles",0,&read_elements_event);
3421 PetscLogEventBegin(read_elements_event,0,0,0,0);
3422 std::ifstream vertexFile(vertexFileName.c_str());
3423 if (!vertexFile.good())
3425 std::cerr<<
"cannot open triangle node file "
3426 <<vertexFileName<<std::endl;
3430 int hasVertexMarkers(0),hasVertexAttributes(0),nSpace(2),nNodes_global;
3432 vertexFile >>
eatcomments >> nNodes_global >> nSpace >> hasVertexAttributes >> hasVertexMarkers >>
eatline ;
3433 assert(nNodes_global > 0);
3434 assert(nSpace == 2);
3439 if (hasVertexAttributes > 0)
3441 std::cerr<<
"WARNING triangle nodes hasAttributes= "<<hasVertexAttributes
3442 <<
" > 0 will treat first value as integer id for boundary!!"<<std::endl;
3443 hasVertexMarkers = 1;
3449 valarray<int> nodeOffsets_old(size+1);
3450 nodeOffsets_old[0] = 0;
3451 for (
int sdN=0; sdN < size; sdN++)
3453 nodeOffsets_old[sdN+1] = nodeOffsets_old[sdN] +
3454 int(nNodes_global)/size + (int(nNodes_global)%size > sdN);
3456 int nNodes_subdomain_old = nodeOffsets_old[rank+1] - nodeOffsets_old[rank];
3463 std::ifstream elementFile(elementFileName.c_str());
3464 if (!elementFile.good())
3466 std::cerr<<
"cannot open triangle element file "
3467 <<elementFileName<<std::endl;
3472 int nNodesPerSimplex(simplexDim),hasElementMarkers = 0,nElements_global;
3473 elementFile >>
eatcomments >> nElements_global >> nNodesPerSimplex >> hasElementMarkers >>
eatline;
3474 assert(nElements_global > 0);
3475 assert(nNodesPerSimplex == simplexDim);
3479 map<int,vector<int> > elements_old;
3480 for (
int ie = 0; ie < nElements_global; ie++)
3485 assert(0 <= ne && ne < nElements_global && elementFile.good());
3486 for (
int iv = 0; iv < simplexDim; iv++)
3490 assert(0 <= nv && nv < nNodes_global);
3491 element_nodes_old[iv] = nv;
3494 for (
int iv = 0; iv < simplexDim; iv++)
3497 int nN_star = element_nodes_old[iv];
3498 bool inSubdomain=
false;
3499 if (nN_star >= nodeOffsets_old[rank] && nN_star < nodeOffsets_old[rank+1])
3503 for (
int jv = 0; jv < simplexDim; jv++)
3507 int nN_star_subdomain = nN_star-nodeOffsets_old[rank];
3508 nodeStar[nN_star_subdomain].insert(element_nodes_old[jv]);
3513 elements_old[ie] = element_nodes_old;
3517 elementFile.close();
3518 PetscLogEventEnd(read_elements_event,0,0,0,0);
3519 int repartition_nodes_event;
3520 PetscLogEventRegister(
"Repart nodes",0,&repartition_nodes_event);
3521 PetscLogEventBegin(repartition_nodes_event,0,0,0,0);
3524 valarray<int> nodeStarOffsets(nNodes_subdomain_old+1);
3525 nodeStarOffsets[0] = 0;
3526 for (
int nN=1;nN<nNodes_subdomain_old+1;nN++)
3527 nodeStarOffsets[nN] = nodeStarOffsets[nN-1] + nodeStar[nN-1].size();
3528 valarray<int> nodeStarArray(nodeStarOffsets[nNodes_subdomain_old]);
3529 for (
int nN=0,offset=0;nN<nNodes_subdomain_old;nN++)
3530 for (set<int>::iterator nN_star=nodeStar[nN].begin();nN_star!=nodeStar[nN].end();nN_star++,offset++)
3531 nodeStarArray[offset] = *nN_star;
3533 int max_nNodeNeighbors_node=0;
3534 for (
int nN=0;nN<nNodes_subdomain_old;nN++)
3535 max_nNodeNeighbors_node=
max(max_nNodeNeighbors_node,nodeStarOffsets[nN+1]-nodeStarOffsets[nN]);
3537 PetscBool isInitialized;
3538 PetscInitialized(&isInitialized);
3539 PetscInt *nodeNeighborsOffsets_subdomain,*nodeNeighbors_subdomain,*weights_subdomain,*vertex_weights_subdomain;
3540 PetscReal *partition_weights;
3541 PetscMalloc(
sizeof(PetscInt)*(nNodes_subdomain_old+1),&nodeNeighborsOffsets_subdomain);
3542 PetscMalloc(
sizeof(PetscInt)*(nNodes_subdomain_old*max_nNodeNeighbors_node),&nodeNeighbors_subdomain);
3543 PetscMalloc(
sizeof(PetscInt)*(nNodes_subdomain_old*max_nNodeNeighbors_node),&weights_subdomain);
3544 PetscMalloc(
sizeof(PetscInt)*(nNodes_subdomain_old),&vertex_weights_subdomain);
3545 PetscMalloc(
sizeof(PetscReal)*(size),&partition_weights);
3546 for (
int sd=0;sd<size;sd++)
3547 partition_weights[sd] = 1.0/
double(size);
3548 nodeNeighborsOffsets_subdomain[0] = 0;
3550 for (
int nN = 0,offset=0; nN < nNodes_subdomain_old; nN++)
3552 for (
int offset_subdomain = nodeStarOffsets[nN];
3553 offset_subdomain < nodeStarOffsets[nN+1];
3556 nodeNeighbors_subdomain[offset++] = nodeStarArray[offset_subdomain];
3558 nodeNeighborsOffsets_subdomain[nN+1]=offset;
3559 sort(&nodeNeighbors_subdomain[nodeNeighborsOffsets_subdomain[nN]],&nodeNeighbors_subdomain[nodeNeighborsOffsets_subdomain[nN+1]]);
3561 int weight= (nodeNeighborsOffsets_subdomain[nN+1] - nodeNeighborsOffsets_subdomain[nN]);
3562 vertex_weights_subdomain[nN] = weight;
3563 for (
int k=nodeNeighborsOffsets_subdomain[nN];k<nodeNeighborsOffsets_subdomain[nN+1];k++)
3564 weights_subdomain[k] = weight;
3570 int nNodes_subdomain_max=0;
3571 MPI_Allreduce(&nNodes_subdomain_old,
3572 &nNodes_subdomain_max,
3576 PROTEUS_COMM_WORLD);
3578 std::cout<<
"Max nNodes_subdomain "<<nNodes_subdomain_max<<
" nNodes_global "<<nNodes_global<<std::endl;
3579 ierr = MatCreateMPIAdj(PROTEUS_COMM_WORLD,
3580 nNodes_subdomain_old,
3582 nodeNeighborsOffsets_subdomain,
3583 nodeNeighbors_subdomain,
3585 &petscAdjacency);CHKERRABORT(PROTEUS_COMM_WORLD, ierr);
3587 const double max_rss_gb(10.0);
3588 ierr =
enforceMemoryLimit(PROTEUS_COMM_WORLD, rank, max_rss_gb,
"Done allocating MPIAdj");CHKERRABORT(PROTEUS_COMM_WORLD, ierr);
3589 MatPartitioning petscPartition;
3590 ierr = MatPartitioningCreate(PROTEUS_COMM_WORLD,&petscPartition);CHKERRABORT(PROTEUS_COMM_WORLD, ierr);
3591 ierr = MatPartitioningSetAdjacency(petscPartition,petscAdjacency);CHKERRABORT(PROTEUS_COMM_WORLD, ierr);
3592 ierr = MatPartitioningSetFromOptions(petscPartition);CHKERRABORT(PROTEUS_COMM_WORLD, ierr);
3593 ierr = MatPartitioningSetVertexWeights(petscPartition,vertex_weights_subdomain);CHKERRABORT(PROTEUS_COMM_WORLD, ierr);
3594 ierr = MatPartitioningSetPartitionWeights(petscPartition,partition_weights);CHKERRABORT(PROTEUS_COMM_WORLD, ierr);
3596 IS nodePartitioningIS_new;
3597 ierr = MatPartitioningApply(petscPartition,&nodePartitioningIS_new);CHKERRABORT(PROTEUS_COMM_WORLD, ierr);
3598 ierr = MatPartitioningDestroy(&petscPartition);CHKERRABORT(PROTEUS_COMM_WORLD, ierr);
3602 valarray<int> nNodes_subdomain_new(size);
3603 ISPartitioningCount(nodePartitioningIS_new,size,&nNodes_subdomain_new[0]);
3606 valarray<int> nodeOffsets_new(size+1);
3607 nodeOffsets_new[0] = 0;
3608 for (
int sdN = 0; sdN < size; sdN++)
3610 nodeOffsets_new[sdN+1] = nodeOffsets_new[sdN] + nNodes_subdomain_new[sdN];
3614 IS nodeNumberingIS_subdomain_old2new;
3615 ISPartitioningToNumbering(nodePartitioningIS_new,&nodeNumberingIS_subdomain_old2new);
3622 hid_t plist_id = H5Pcreate(H5P_FILE_ACCESS);
3623#ifdef H5_HAVE_PARALLEL
3624 H5Pset_fapl_mpio(plist_id, PROTEUS_COMM_WORLD, MPI_INFO_NULL);
3630 const char* H5FILE_NAME(
"mappings.h5");
3631 hid_t file_id = H5Fcreate(H5FILE_NAME, H5F_ACC_TRUNC, H5P_DEFAULT, plist_id);
3639 dimsf[0] = nNodes_global;
3641 hid_t filespace = H5Screate_simple(
RANK, dimsf, NULL);
3646 hid_t dset_id = H5Dcreate(file_id,
"nodeNumbering_old2new", H5T_NATIVE_INT, filespace,
3647 H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
3648 H5Sclose(filespace);
3656 count[0] = nNodes_subdomain_old;
3657 offset[0] = nodeOffsets_old[rank];
3658 hid_t memspace = H5Screate_simple(
RANK, count, NULL);
3663 filespace = H5Dget_space(dset_id);
3664 H5Sselect_hyperslab(filespace, H5S_SELECT_SET, offset, NULL, count, NULL);
3673 const PetscInt* data;
3674 ISGetIndices(nodeNumberingIS_subdomain_old2new, &data);
3679 plist_id = H5Pcreate(H5P_DATASET_XFER);
3680#ifdef H5_HAVE_PARALLEL
3681 H5Pset_dxpl_mpio(plist_id, H5FD_MPIO_COLLECTIVE);
3684 herr_t status = H5Dwrite(dset_id, H5T_NATIVE_INT, memspace, filespace,
3687 ISRestoreIndices(nodeNumberingIS_subdomain_old2new, &data);
3698 IS nodeNumberingIS_global_old2new;
3699 ISAllGather(nodeNumberingIS_subdomain_old2new,&nodeNumberingIS_global_old2new);
3700 const PetscInt * nodeNumbering_global_old2new;
3701 ISGetIndices(nodeNumberingIS_global_old2new,&nodeNumbering_global_old2new);
3709 int dset_data[nNodes_global];
3715 dataset_id = H5Dopen2(file_id,
"/nodeNumbering_old2new", H5P_DEFAULT);
3717 status = H5Dread(dataset_id, H5T_NATIVE_INT, H5S_ALL, H5S_ALL, H5P_DEFAULT,
3721 status = H5Dclose(dataset_id);
3723 for (
int i=0;i<nNodes_global;i++)
3724 assert(nodeNumbering_global_old2new[i] == dset_data[i]);
3725 std::cout<<
"==================out of core old2new is correct!===================="<<std::endl;
3737 PetscLogEventEnd(repartition_nodes_event,0,0,0,0);
3738 int receive_element_mask_event;
3739 PetscLogEventRegister(
"Recv. ele mask",0,&receive_element_mask_event);
3740 PetscLogEventBegin(receive_element_mask_event,0,0,0,0);
3745 PetscLogEventEnd(receive_element_mask_event,0,0,0,0);
3747 int build_subdomains_reread_elements_event;
3748 PetscLogEventRegister(
"Reread eles",0,&build_subdomains_reread_elements_event);
3749 PetscLogEventBegin(build_subdomains_reread_elements_event,0,0,0,0);
3756 std::ifstream elementFile2(elementFileName.c_str());
3758 if (!elementFile2.good())
3760 std::cerr<<
"cannot open triangle elements file"
3761 <<elementFileName<<std::endl;
3765 elementFile2 >>
eatcomments >> nElements_global >> nNodesPerSimplex >> hasElementMarkers >>
eatline;
3766 assert(nElements_global > 0);
3767 assert(nNodesPerSimplex == simplexDim);
3768 set<int> elements_subdomain_owned;
3770 int element_nodes_new_array[3];
3773 map<int,vector<int> > elementNodesArrayMap;
3774 map<int,long int> elementMaterialTypesMap;
3776 map<NodeTuple<2>,set<pair<int,int> > > edgeElementsMap;
3779 for (
int ie = 0; ie < nElements_global; ie++)
3781 int ne, nv, elementId(0);
3782 long double elementId_double;
3785 assert(0 <= ne && ne < nElements_global && elementFile.good());
3786 for (
int iv = 0; iv < simplexDim; iv++)
3788 elementFile2 >> nv ;
3790 assert(0 <= nv && nv < nNodes_global);
3791 element_nodes_old[iv] = nv;
3792 element_nodes_new[iv] = nodeNumbering_global_old2new[nv];
3793 element_nodes_new_array[iv] = element_nodes_new[iv];
3796 for (
int iv = 0; iv < simplexDim; iv++)
3798 int nN_star_new = element_nodes_new[iv];
3799 bool inSubdomain=
false;
3800 if (nN_star_new >= nodeOffsets_new[rank] && nN_star_new < nodeOffsets_new[rank+1])
3804 for (
int ebN=0;ebN < simplexDim ; ebN++)
3806 int nodes[simplexDim-1] = { element_nodes_new[(ebN+1) % simplexDim],
3807 element_nodes_new[(ebN+2) % simplexDim]};
3808 NodeTuple<simplexDim-1> nodeTuple(nodes);
3809 if(elementBoundaryElementsMap.find(nodeTuple) != elementBoundaryElementsMap.end())
3811 if (elementBoundaryElementsMap[nodeTuple].right == -1 && ne != elementBoundaryElementsMap[nodeTuple].left)
3813 elementBoundaryElementsMap[nodeTuple].right=ne;
3814 elementBoundaryElementsMap[nodeTuple].right_ebN_element=ebN;
3823 for (
int nNL=0,edN=0;nNL < simplexDim ; nNL++)
3824 for(
int nNR=nNL+1;nNR < simplexDim;nNR++,edN++)
3826 int nodes[2] = { element_nodes_new[nNL],
3827 element_nodes_new[nNR]};
3829 edgeElementsMap[nodeTuple].insert(pair<int,int>(ne,edN));
3832 int nN_star_new_subdomain = nN_star_new - nodeOffsets_new[rank];
3833 nodeElementsStar[nN_star_new_subdomain].insert(ne);
3834 for (
int jv = 0; jv < simplexDim; jv++)
3838 int nN_point_new = element_nodes_new[jv];
3839 nodeStarNew[nN_star_new_subdomain].insert(nN_point_new);
3845 elementNodesArrayMap[ne] = element_nodes_new;
3848 if (elementNodesArrayMap.find(ne) != elementNodesArrayMap.end())
3850 if (nodeTuple.
nodes[1] >= nodeOffsets_new[rank] && nodeTuple.
nodes[1] < nodeOffsets_new[rank+1])
3851 elements_subdomain_owned.insert(ne);
3852 if (hasElementMarkers > 0)
3854 elementFile2 >> elementId_double;
3855 elementId =
static_cast<long int>(elementId_double);
3856 elementMaterialTypesMap[ne] = elementId;
3861 elementFile2.close();
3862 int nElements_owned_subdomain(elements_subdomain_owned.size()),
3863 nElements_owned_new=0;
3864 MPI_Allreduce(&nElements_owned_subdomain,&nElements_owned_new,1,MPI_INT,MPI_SUM,PROTEUS_COMM_WORLD);
3865 assert(nElements_owned_new == nElements_global);
3866 PetscLogEventEnd(build_subdomains_reread_elements_event,0,0,0,0);
3867 int build_subdomains_send_marked_elements_event;
3868 PetscLogEventRegister(
"Mark/send eles",0,&build_subdomains_send_marked_elements_event);
3869 PetscLogEventBegin(build_subdomains_send_marked_elements_event,0,0,0,0);
3875 valarray<int> nodeElementOffsets(nNodes_subdomain_new[rank]+1);
3876 nodeElementOffsets[0] = 0;
3877 for (
int nN = 0; nN < nNodes_subdomain_new[rank]; nN++)
3878 nodeElementOffsets[nN+1] = nodeElementOffsets[nN]+nodeElementsStar[nN].size();
3879 valarray<int> nodeElementsArray(nodeElementOffsets[nNodes_subdomain_new[rank]]);
3880 for (
int nN=0,offset=0; nN < nNodes_subdomain_new[rank]; nN++)
3882 for (set<int>::iterator eN_star = nodeElementsStar[nN].begin(); eN_star != nodeElementsStar[nN].end();
3885 nodeElementsArray[offset] = *eN_star;
3889 valarray<int> nodeStarOffsetsNew(nNodes_subdomain_new[rank]+1);
3890 nodeStarOffsetsNew[0] = 0;
3891 for (
int nN=1;nN<nNodes_subdomain_new[rank]+1;nN++)
3892 nodeStarOffsetsNew[nN] = nodeStarOffsetsNew[nN-1] + nodeStarNew[nN-1].size();
3893 valarray<int> nodeStarArrayNew(nodeStarOffsetsNew[nNodes_subdomain_new[rank]]);
3894 for (
int nN=0,offset=0;nN<nNodes_subdomain_new[rank];nN++)
3896 for (set<int>::iterator nN_star=nodeStarNew[nN].begin();nN_star!=nodeStarNew[nN].end();nN_star++,offset++)
3898 nodeStarArrayNew[offset] = *nN_star;
3901 PetscLogEventEnd(build_subdomains_send_marked_elements_event,0,0,0,0);
3902 int build_subdomains_global_numbering_elements_event;
3903 PetscLogEventRegister(
"Global ele nmbr",0,&build_subdomains_global_numbering_elements_event);
3904 PetscLogEventBegin(build_subdomains_global_numbering_elements_event,0,0,0,0);
3908 valarray<int> nElements_subdomain_new(size),
3909 elementOffsets_new(size+1);
3910 for (
int sdN = 0; sdN < size; sdN++)
3913 nElements_subdomain_new[sdN] = int(elements_subdomain_owned.size());
3915 nElements_subdomain_new[sdN] = 0;
3917 valarray<int> nElements_subdomain_new_send = nElements_subdomain_new;
3918 MPI_Allreduce(&nElements_subdomain_new_send[0],&nElements_subdomain_new[0],size,MPI_INT,MPI_SUM,PROTEUS_COMM_WORLD);
3920 elementOffsets_new[0] = 0;
3921 for (
int sdN = 0; sdN < size; sdN++)
3922 elementOffsets_new[sdN+1] = elementOffsets_new[sdN] + nElements_subdomain_new[sdN];
3924 valarray<int> elementNumbering_subdomain_new2old(elements_subdomain_owned.size());
3925 set<int>::iterator eN_ownedp = elements_subdomain_owned.begin();
3926 for (
int eN = 0; eN < int(elements_subdomain_owned.size()); eN++,eN_ownedp++)
3928 elementNumbering_subdomain_new2old[eN] = *eN_ownedp;
3931 IS elementNumberingIS_subdomain_new2old;
3932 ISCreateGeneral(PROTEUS_COMM_WORLD,elements_subdomain_owned.size(),&elementNumbering_subdomain_new2old[0],PETSC_COPY_VALUES,
3933 &elementNumberingIS_subdomain_new2old);
3934 IS elementNumberingIS_global_new2old;
3935 ISAllGather(elementNumberingIS_subdomain_new2old,&elementNumberingIS_global_new2old);
3937 const PetscInt *elementNumbering_global_new2old;
3938 ISGetIndices(elementNumberingIS_global_new2old,&elementNumbering_global_new2old);
3940 valarray<int> elementNumbering_global_old2new(nElements_global);
3941 for (
int eN = 0; eN < nElements_global; eN++)
3943 elementNumbering_global_old2new[elementNumbering_global_new2old[eN]] = eN;
3945 PetscLogEventEnd(build_subdomains_global_numbering_elements_event,0,0,0,0);
3947 int build_subdomains_faces_event;
3948 PetscLogEventRegister(
"Subd faces",0,&build_subdomains_faces_event);
3949 PetscLogEventBegin(build_subdomains_faces_event,0,0,0,0);
3961 std::ifstream elementBoundaryFile(elementBoundaryFileName.c_str());
3963 if (!elementBoundaryFile.good())
3965 std::cerr<<
"cannot open triangle edge file "
3966 <<elementBoundaryFileName<<std::endl;
3971 bool hasElementBoundaryMarkers =
false;
3972 int nElementBoundaries_global;
3973 int ihasElementBoundaryMarkers(0);
3975 elementBoundaryFile >>
eatcomments >> nElementBoundaries_global >>ihasElementBoundaryMarkers >>
eatline ;
3976 assert(nElementBoundaries_global > 0);
3977 if (ihasElementBoundaryMarkers > 0)
3979 hasElementBoundaryMarkers =
true;
3983 set<int> elementBoundaries_subdomain_owned;
3984 vector<set<int> > nodeElementBoundariesStar(nNodes_subdomain_new[rank]);
3985 map<int,int> elementBoundaryMaterialTypesMap;
3986 map<int,vector<int> > elementBoundariesMap;
3987 set<int> supportedElementBoundaries;
3988 for (
int ieb = 0; ieb < nElementBoundaries_global; ieb++)
3990 int neb,nn0,nn1;
int ebId(0);
3991 elementBoundaryFile >>
eatcomments >> neb >> nn0 >> nn1;
3992 if (ihasElementBoundaryMarkers > 0)
3993 elementBoundaryFile >> ebId;
3997 assert(0 <= neb && neb < nElementBoundaries_global && elementBoundaryFile.good());
4000 int nn0_new = nodeNumbering_global_old2new[nn0];
4001 if (nn0_new >= nodeOffsets_new[rank] && nn0_new < nodeOffsets_new[rank+1])
4003 nodeElementBoundariesStar[nn0_new-nodeOffsets_new[rank]].insert(neb);
4004 supportedElementBoundaries.insert(neb);
4006 int nn1_new = nodeNumbering_global_old2new[nn1];
4007 if (nn1_new >= nodeOffsets_new[rank] && nn1_new < nodeOffsets_new[rank+1])
4009 nodeElementBoundariesStar[nn1_new-nodeOffsets_new[rank]].insert(neb);
4010 supportedElementBoundaries.insert(neb);
4012 int nodes[2] = {nn0_new,nn1_new};
4014 elementBoundaryFile >>
eatline;
4015 if (elementBoundaryElementsMap.find(nodeTuple) != elementBoundaryElementsMap.end())
4017 if (nodeTuple.
nodes[1] >= nodeOffsets_new[rank] && nodeTuple.
nodes[1] < nodeOffsets_new[rank+1])
4018 elementBoundaries_subdomain_owned.insert(neb);
4019 if (ihasElementBoundaryMarkers > 0)
4020 elementBoundaryMaterialTypesMap[neb]=ebId;
4021 int eN_left = elementNumbering_global_old2new[elementBoundaryElementsMap[nodeTuple].left];
4022 if (elementBoundariesMap.find(eN_left) != elementBoundariesMap.end())
4024 elementBoundariesMap[eN_left][elementBoundaryElementsMap[nodeTuple].left_ebN_element] = neb;
4030 elementBoundariesMap[eN_left] = elementBoundaries_element;
4032 elementBoundariesMap[eN_left][elementBoundaryElementsMap[nodeTuple].left_ebN_element] = neb;
4034 if (elementBoundaryElementsMap[nodeTuple].right >= 0)
4036 int eN_right = elementNumbering_global_old2new[elementBoundaryElementsMap[nodeTuple].right];
4037 if (elementBoundariesMap.find(eN_right) != elementBoundariesMap.end())
4039 elementBoundariesMap[eN_right][elementBoundaryElementsMap[nodeTuple].right_ebN_element] = neb;
4045 elementBoundariesMap[eN_right] = elementBoundaries_element;
4047 elementBoundariesMap[eN_right][elementBoundaryElementsMap[nodeTuple].right_ebN_element] = neb;
4052 elementBoundaryElementsMap.clear();
4054 elementBoundaryFile.close();
4056 int nElementBoundaries_owned_subdomain=elementBoundaries_subdomain_owned.size(),
4057 nElementBoundaries_owned_new=0;
4058 MPI_Allreduce(&nElementBoundaries_owned_subdomain,&nElementBoundaries_owned_new,1,MPI_INT,MPI_SUM,PROTEUS_COMM_WORLD);
4059 assert(nElementBoundaries_owned_new == nElementBoundaries_global);
4062 for (map<
int,
vector<int> >::iterator elementBoundariesp=elementBoundariesMap.begin();
4063 elementBoundariesp!=elementBoundariesMap.end();
4064 elementBoundariesp++)
4067 for (
int iv=0;iv<simplexDim;iv++)
4070 int nN_global = elementNodesArrayMap[elementNumbering_global_new2old[elementBoundariesp->first]][iv];
4071 if (nN_global >= nodeOffsets_new[rank] && nN_global < nodeOffsets_new[rank+1])
4074 for(
int eb=0;eb<simplexDim;eb++)
4076 nodeElementBoundariesStar[nN_global-nodeOffsets_new[rank]].insert(elementBoundariesp->second[eb]);
4082 valarray<int> nodeElementBoundaryOffsets(nNodes_subdomain_new[rank]+1);
4083 nodeElementBoundaryOffsets[0] = 0;
4084 for (
int nN = 0; nN < nNodes_subdomain_new[rank]; nN++)
4085 nodeElementBoundaryOffsets[nN+1] = nodeElementBoundaryOffsets[nN]+nodeElementBoundariesStar[nN].size();
4086 valarray<int> nodeElementBoundariesArray(nodeElementBoundaryOffsets[nNodes_subdomain_new[rank]]);
4087 for (
int nN=0,offset=0; nN < nNodes_subdomain_new[rank]; nN++)
4089 for (set<int>::iterator ebN_star = nodeElementBoundariesStar[nN].begin(); ebN_star != nodeElementBoundariesStar[nN].end();
4090 ebN_star++,offset++)
4092 nodeElementBoundariesArray[offset] = *ebN_star;
4097 valarray<int> nElementBoundaries_subdomain_new(size),
4098 elementBoundaryOffsets_new(size+1);
4099 for (
int sdN=0;sdN<size;sdN++)
4101 nElementBoundaries_subdomain_new[sdN] = elementBoundaries_subdomain_owned.size();
4103 nElementBoundaries_subdomain_new[sdN] = 0;
4104 valarray<int> nElementBoundaries_subdomain_new_send=nElementBoundaries_subdomain_new;
4105 MPI_Allreduce(&nElementBoundaries_subdomain_new_send[0],&nElementBoundaries_subdomain_new[0],size,MPI_INT,MPI_SUM,PROTEUS_COMM_WORLD);
4106 elementBoundaryOffsets_new[0] = 0;
4107 for (
int sdN=0;sdN<size;sdN++)
4108 elementBoundaryOffsets_new[sdN+1] = elementBoundaryOffsets_new[sdN]+nElementBoundaries_subdomain_new[sdN];
4114 valarray<int> elementBoundaryNumbering_new2old(elementBoundaries_subdomain_owned.size());
4115 set<int>::iterator ebN_ownedp=elementBoundaries_subdomain_owned.begin();
4116 for (
int ebN=0;ebN<int(elementBoundaries_subdomain_owned.size());ebN++)
4118 elementBoundaryNumbering_new2old[ebN] = *ebN_ownedp++;
4120 IS elementBoundaryNumberingIS_subdomain_new2old;
4121 ISCreateGeneral(PROTEUS_COMM_WORLD,elementBoundaries_subdomain_owned.size(),&elementBoundaryNumbering_new2old[0],PETSC_COPY_VALUES,&elementBoundaryNumberingIS_subdomain_new2old);
4122 IS elementBoundaryNumberingIS_global_new2old;
4123 ISAllGather(elementBoundaryNumberingIS_subdomain_new2old,&elementBoundaryNumberingIS_global_new2old);
4124 const PetscInt *elementBoundaryNumbering_global_new2old;
4126 ISGetIndices(elementBoundaryNumberingIS_global_new2old,&elementBoundaryNumbering_global_new2old);
4130 elementBoundaryNumbering_global_old2new[elementBoundaryNumbering_global_new2old[ebN]] = ebN;
4132 ISRestoreIndices(elementBoundaryNumberingIS_global_new2old,&elementBoundaryNumbering_global_new2old);
4133 ISDestroy(&elementBoundaryNumberingIS_subdomain_new2old);
4134 ISDestroy(&elementBoundaryNumberingIS_global_new2old);
4135 PetscLogEventEnd(build_subdomains_faces_event,0,0,0,0);
4137 int build_subdomains_edges_event;
4138 PetscLogEventRegister(
"Subd edges",0,&build_subdomains_edges_event);
4139 PetscLogEventBegin(build_subdomains_edges_event,0,0,0,0);
4144 std::ifstream edgeFile(edgeFileName.c_str());
4146 if (!edgeFile.good())
4148 std::cerr<<
"cannot open Triangle edge file"
4149 <<edgeFileName<<std::endl;
4154 bool hasEdgeMarkers =
false;
4156 int ihasEdgeMarkers(0);
4158 assert(nEdges_global > 0);
4159 if (ihasEdgeMarkers > 0)
4161 hasEdgeMarkers =
true;
4165 set<int> edges_subdomain_owned;
4167 map<int,int> edgeMaterialTypesMap;
4168 map<int,vector<int> > elementEdgesMap;
4169 map<int,pair<int,int> > edgeNodesMap;
4170 set<int> supportedEdges;
4171 for (
int ied = 0; ied < nEdges_global; ied++)
4173 int ned,nn0,nn1;
int edId(0);
4175 if (ihasEdgeMarkers > 0)
4180 assert(0 <= ned && ned < nEdges_global && edgeFile.good());
4181 int nn0_new = nodeNumbering_global_old2new[nn0];
4182 if (nn0_new >= nodeOffsets_new[rank] && nn0_new < nodeOffsets_new[rank+1])
4184 nodeEdgesStar.at(nn0_new-nodeOffsets_new[rank]).insert(ned);
4185 supportedEdges.insert(ned);
4187 int nn1_new = nodeNumbering_global_old2new[nn1];
4188 if (nn1_new >= nodeOffsets_new[rank] && nn1_new < nodeOffsets_new[rank+1])
4190 nodeEdgesStar.at(nn1_new-nodeOffsets_new[rank]).insert(ned);
4191 supportedEdges.insert(ned);
4193 int nodes[2] = {nn0_new,nn1_new};
4196 if (edgeElementsMap.find(nodeTuple) != edgeElementsMap.end())
4198 if (nodeTuple.
nodes[0] >= nodeOffsets_new[rank] && nodeTuple.
nodes[0] < nodeOffsets_new[rank+1])
4199 edges_subdomain_owned.insert(ned);
4201 edgeNodesMap[ned].first = nodeTuple.
nodes[0];
4202 edgeNodesMap[ned].second = nodeTuple.
nodes[1];
4203 if (ihasEdgeMarkers > 0)
4204 edgeMaterialTypesMap[ned]=edId;
4205 for (set<pair<int,int> >::iterator elementp=edgeElementsMap[nodeTuple].begin();
4206 elementp != edgeElementsMap[nodeTuple].end();
4209 int eN = elementNumbering_global_old2new[elementp->first];
4210 if (elementEdgesMap.find(eN) != elementEdgesMap.end())
4212 elementEdgesMap[eN][elementp->second] = ned;
4216 std::vector<int>
init(3,-1);
4217 elementEdgesMap[eN] =
init;
4218 elementEdgesMap[eN][elementp->second] = ned;
4223 edgeElementsMap.clear();
4225 int nEdges_owned_subdomain=edges_subdomain_owned.size(),
4227 MPI_Allreduce(&nEdges_owned_subdomain,&nEdges_owned_new,1,MPI_INT,MPI_SUM,PROTEUS_COMM_WORLD);
4228 assert(nEdges_owned_new == nEdges_global);
4235 for (map<
int,
vector<int> >::iterator edgesp=elementEdgesMap.begin();
4236 edgesp!=elementEdgesMap.end();
4240 for (
int iv=0;iv<simplexDim;iv++)
4243 int nN_global = elementNodesArrayMap[elementNumbering_global_new2old[edgesp->first]][iv];
4244 if (nN_global >= nodeOffsets_new[rank] && nN_global < nodeOffsets_new[rank+1])
4247 for(
int ed=0;ed<3;ed++)
4249 nodeEdgesStar.at(nN_global-nodeOffsets_new[rank]).insert(edgesp->second[ed]);
4255 valarray<int> nodeEdgeOffsets(nNodes_subdomain_new[rank]+1);
4256 nodeEdgeOffsets[0] = 0;
4257 for (
int nN = 0; nN < nNodes_subdomain_new[rank]; nN++)
4258 nodeEdgeOffsets[nN+1] = nodeEdgeOffsets[nN]+nodeEdgesStar.at(nN).size();
4259 valarray<int> nodeEdgesArray(nodeEdgeOffsets[nNodes_subdomain_new[rank]]);
4260 for (
int nN=0,offset=0; nN < nNodes_subdomain_new[rank]; nN++)
4262 for (set<int>::iterator edN_star = nodeEdgesStar.at(nN).begin();
4263 edN_star != nodeEdgesStar.at(nN).end();
4264 edN_star++,offset++)
4266 nodeEdgesArray[offset] = *edN_star;
4271 valarray<int> nEdges_subdomain_new(size),
4272 edgeOffsets_new(size+1);
4273 for (
int sdN=0;sdN<size;sdN++)
4275 nEdges_subdomain_new[sdN] = edges_subdomain_owned.size();
4277 nEdges_subdomain_new[sdN] = 0;
4278 valarray<int> nEdges_subdomain_new_send=nEdges_subdomain_new;
4279 MPI_Allreduce(&nEdges_subdomain_new_send[0],&nEdges_subdomain_new[0],size,MPI_INT,MPI_SUM,PROTEUS_COMM_WORLD);
4283 edgeOffsets_new[0] = 0;
4284 for (
int sdN=0;sdN<size;sdN++)
4285 edgeOffsets_new[sdN+1] = edgeOffsets_new[sdN]+nEdges_subdomain_new[sdN];
4291 valarray<int> edgeNumbering_new2old(edges_subdomain_owned.size());
4292 set<int>::iterator edN_ownedp=edges_subdomain_owned.begin();
4293 for (
int edN=0;edN<int(edges_subdomain_owned.size());edN++,edN_ownedp++)
4295 edgeNumbering_new2old[edN] = *edN_ownedp;
4297 IS edgeNumberingIS_subdomain_new2old;
4298 ISCreateGeneral(PROTEUS_COMM_WORLD,edges_subdomain_owned.size(),&edgeNumbering_new2old[0],PETSC_COPY_VALUES,&edgeNumberingIS_subdomain_new2old);
4299 IS edgeNumberingIS_global_new2old;
4300 ISAllGather(edgeNumberingIS_subdomain_new2old,&edgeNumberingIS_global_new2old);
4301 const PetscInt *edgeNumbering_global_new2old;
4302 valarray<int> edgeNumbering_global_old2new(newMesh.
nEdges_global);
4303 ISGetIndices(edgeNumberingIS_global_new2old,&edgeNumbering_global_new2old);
4307 edgeNumbering_global_old2new[edgeNumbering_global_new2old[edN]] = edN;
4309 ISRestoreIndices(edgeNumberingIS_global_new2old,&edgeNumbering_global_new2old);
4310 ISDestroy(&edgeNumberingIS_subdomain_new2old);
4311 ISDestroy(&edgeNumberingIS_global_new2old);
4317 set<int> elements_overlap,nodes_overlap,elementBoundaries_overlap,edges_overlap;
4318 for (
int nN = 0; nN < nNodes_subdomain_new[rank]; nN++)
4321 for (
int offset = nodeStarOffsetsNew[nN];offset<nodeStarOffsetsNew[nN+1];offset++)
4323 int nN_point_global = nodeStarArrayNew[offset];
4324 bool offproc = nN_point_global < nodeOffsets_new[rank] || nN_point_global >= nodeOffsets_new[rank+1];
4326 nodes_overlap.insert(nN_point_global);
4329 for (
int eN_star_offset = nodeElementOffsets[nN];
4330 eN_star_offset < nodeElementOffsets[nN+1]; eN_star_offset++)
4332 int eN_star_old = nodeElementsArray[eN_star_offset];
4333 int eN_star_new = elementNumbering_global_old2new[eN_star_old];
4334 bool offproc = eN_star_new >= elementOffsets_new[rank+1] || eN_star_new < elementOffsets_new[rank];
4336 elements_overlap.insert(eN_star_new);
4339 for (
int ebN_star_offset = nodeElementBoundaryOffsets[nN];
4340 ebN_star_offset < nodeElementBoundaryOffsets[nN+1]; ebN_star_offset++)
4342 int ebN_star_old = nodeElementBoundariesArray[ebN_star_offset];
4343 int ebN_star_new = elementBoundaryNumbering_global_old2new[ebN_star_old];
4344 bool offproc = ebN_star_new >= elementBoundaryOffsets_new[rank+1] || ebN_star_new < elementBoundaryOffsets_new[rank];
4346 elementBoundaries_overlap.insert(ebN_star_new);
4349 for (
int edN_star_offset = nodeEdgeOffsets[nN];
4350 edN_star_offset < nodeEdgeOffsets[nN+1]; edN_star_offset++)
4352 int edN_star_old = nodeEdgesArray[edN_star_offset];
4353 int edN_star_new = edgeNumbering_global_old2new[edN_star_old];
4354 bool offproc = edN_star_new >= edgeOffsets_new[rank+1] || edN_star_new < edgeOffsets_new[rank];
4356 edges_overlap.insert(edN_star_new);
4359 elementNumbering_global_old2new.resize(0);
4360 MPI_Barrier(PROTEUS_COMM_WORLD);
4361 assert(edges_overlap.size() + nEdges_subdomain_new[rank] == edgeNodesMap.size());
4365 int nN_subdomain = nNodes_subdomain_new[rank];
4366 map<int,int> nodes_overlap_global2subdomainMap;
4367 for (set<int>::iterator nN_globalp=nodes_overlap.begin();nN_globalp != nodes_overlap.end(); nN_globalp++,nN_subdomain++)
4368 nodes_overlap_global2subdomainMap[*nN_globalp] = nN_subdomain;
4387 PetscLogEventEnd(build_subdomains_edges_event,0,0,0,0);
4388 int build_subdomains_renumber_event;
4389 PetscLogEventRegister(
"Subd's renumber",0,&build_subdomains_renumber_event);
4390 PetscLogEventBegin(build_subdomains_renumber_event,0,0,0,0);
4395 std::cerr<<
"USER WARNING: In order to avoid a segmentation fault, you need to have supplied the 'f' flag to the triangleOptions input."<<std::endl;
4396 std::cerr<<
"USER WARNING: In order to avoid an edge assertion error, you need to have supplied the 'ee' flag to the triangleOptions input."<<std::endl;
4414 map<int,int> nodeNumbering_global2subdomainMap;
4415 map<int,int> elementBoundaryNumbering_global2subdomainMap;
4416 map<int,int> edgeNumbering_global2subdomainMap;
4422 for (
int iv = 0; iv < nNodes_global; iv++)
4424 int nv;
double x,y;
int nodeId(0);
4426 if (hasVertexMarkers > 0)
4427 vertexFile >> nodeId;
4429 assert(0 <= nv && nv < nNodes_global && vertexFile.good());
4430 int nN_global_new = nodeNumbering_global_old2new[nv];
4432 if (nN_global_new >= nodeOffsets_new[rank] && nN_global_new < nodeOffsets_new[rank+1])
4434 int nv_subdomain_new = nN_global_new - nodeOffsets_new[rank];
4435 nodeNumbering_subdomain2global[nv_subdomain_new] = nN_global_new;
4436 nodeNumbering_global2subdomainMap[nN_global_new] = nv_subdomain_new;
4440 if (hasVertexMarkers > 0)
4444 if (nodes_overlap.count(nN_global_new) == 1)
4446 int nv_subdomain_new = nodes_overlap_global2subdomainMap[nN_global_new];
4447 nodeNumbering_subdomain2global[nv_subdomain_new] = nN_global_new;
4448 nodeNumbering_global2subdomainMap[nN_global_new] = nv_subdomain_new;
4452 if (hasVertexMarkers > 0)
4458 ISRestoreIndices(nodeNumberingIS_global_old2new,&nodeNumbering_global_old2new);
4459 ISDestroy(&nodePartitioningIS_new);
4460 ISDestroy(&nodeNumberingIS_subdomain_old2new);
4461 ISDestroy(&nodeNumberingIS_global_old2new);
4472 for (
int eN = 0; eN < nElements_subdomain_new[rank]; eN++)
4474 int eN_global_new = elementOffsets_new[rank] + eN;
4475 int eN_global_old = elementNumbering_global_new2old[eN_global_new];
4476 elementNumbering_subdomain2global[eN] = eN_global_new;
4480 int nN_global_new = elementNodesArrayMap[eN_global_old][nN];
4481 int nN_subdomain = nodeNumbering_global2subdomainMap[nN_global_new];
4488 set<int>::iterator eN_p = elements_overlap.begin();
4489 for (
int eN = nElements_subdomain_new[rank]; eN < nElements_subdomain_new[rank] + int(elements_overlap.size()); eN++,eN_p++)
4491 int eN_global_new = *eN_p;
4492 int eN_global_old = elementNumbering_global_new2old[eN_global_new];
4493 elementNumbering_subdomain2global[eN] = eN_global_new;
4497 int nN_global_new = elementNodesArrayMap[eN_global_old][nN];
4498 int nN_subdomain = nodeNumbering_global2subdomainMap[nN_global_new];
4502 elementNodesArrayMap.clear();
4503 elementMaterialTypesMap.clear();
4504 ISRestoreIndices(elementNumberingIS_global_new2old,&elementNumbering_global_new2old);
4505 ISDestroy(&elementNumberingIS_subdomain_new2old);
4506 ISDestroy(&elementNumberingIS_global_new2old);
4512 for (
int ebN=0; ebN < nElementBoundaries_subdomain_new[rank]; ebN++)
4514 int ebN_global = ebN + elementBoundaryOffsets_new[rank];
4515 elementBoundaryNumbering_subdomain2global[ebN]=ebN_global;
4516 elementBoundaryNumbering_global2subdomainMap[ebN_global] = ebN;
4521 set<int>::iterator ebN_p = elementBoundaries_overlap.begin();
4522 for(
int ebN=nElementBoundaries_subdomain_new[rank];ebN < nElementBoundaries_subdomain_new[rank] + int(elementBoundaries_overlap.size()); ebN++,ebN_p++)
4524 int ebN_global = *ebN_p;
4525 elementBoundaryNumbering_subdomain2global[ebN] = ebN_global;
4526 elementBoundaryNumbering_global2subdomainMap[ebN_global] = ebN;
4535 for (
int eN=0;eN<nElements_subdomain_new[rank];eN++)
4537 int eN_global = eN+elementOffsets_new[rank];
4541 elementBoundaryNumbering_global2subdomainMap[elementBoundaryNumbering_global_old2new[elementBoundariesMap[eN_global][ebN]]];
4547 eN_p = elements_overlap.begin();
4548 for (
int eN = nElements_subdomain_new[rank]; eN < nElements_subdomain_new[rank] + int(elements_overlap.size()); eN++,eN_p++)
4550 int eN_global_new = *eN_p;
4554 elementBoundaryNumbering_global2subdomainMap[elementBoundaryNumbering_global_old2new[elementBoundariesMap[eN_global_new][ebN]]];
4562 for (
int edN=0; edN < nEdges_subdomain_new[rank]; edN++)
4564 int edN_global = edN + edgeOffsets_new[rank];
4565 edgeNumbering_subdomain2global[edN]=edN_global;
4566 edgeNumbering_global2subdomainMap[edN_global] = edN;
4571 set<int>::iterator edN_p = edges_overlap.begin();
4572 for(
int edN=nEdges_subdomain_new[rank];edN < nEdges_subdomain_new[rank] + int(edges_overlap.size()); edN++,edN_p++)
4574 int edN_global = *edN_p;
4575 edgeNumbering_subdomain2global[edN] = edN_global;
4576 edgeNumbering_global2subdomainMap[edN_global] = edN;
4584 for (map<
int,pair<int,int> >::iterator edgep=edgeNodesMap.begin();
4585 edgep!=edgeNodesMap.end();
4588 int edN_global_old = edgep->first;
4589 int edN_global_new = edgeNumbering_global_old2new[edN_global_old];
4590 assert(edgeNumbering_global2subdomainMap.find(edN_global_new) != edgeNumbering_global2subdomainMap.end());
4591 int edN_subdomain = edgeNumbering_global2subdomainMap[edN_global_new];
4595 edgeNumbering_global_old2new.resize(0);
4616 using namespace std;
4621 int> elementBoundaryIds;
4633 if(elementBoundaryElements.find(ebt) != elementBoundaryElements.end())
4635 elementBoundaryElements[ebt].right=eN;
4636 elementBoundaryElements[ebt].right_ebN_element=ebN;
4637 assert(elementBoundaryIds[ebt] == ebN_global);
4641 elementBoundaryElements.insert(elementBoundaryElements.end(),make_pair(ebt,
ElementNeighbors(eN,ebN)));
4642 elementBoundaryIds.insert(elementBoundaryIds.end(),make_pair(ebt,ebN_global));
4652 set<int> interiorElementBoundaries,exteriorElementBoundaries;
4663 eb != elementBoundaryElements.end();
4666 int ebN = elementBoundaryIds[eb->first];
4676 if(eb->second.right != -1)
4678 interiorElementBoundaries.insert(ebN);
4682 exteriorElementBoundaries.insert(ebN);
4684 if (eb->second.right != -1)
4694 for (set<int>::iterator ebN=interiorElementBoundaries.begin();ebN != interiorElementBoundaries.end(); ebN++,ebNI++)
4696 for (set<int>::iterator ebN=exteriorElementBoundaries.begin();ebN != exteriorElementBoundaries.end(); ebN++,ebNE++)
4698 set<NodeTuple<2> > edges;
4723 for (set<int>::iterator nN_star=nodeStar[nN].begin();nN_star!=nodeStar[nN].end();nN_star++,offset++)
4743 for (set<int>::iterator eN_star = nodeElementsStar[nN].begin(); eN_star != nodeElementsStar[nN].end();
4782 if (hasElementBoundaryMarkers)
4785 for (map<int,int>::iterator ebmp = elementBoundaryMaterialTypesMap.begin(); ebmp != elementBoundaryMaterialTypesMap.end();ebmp++)
4787 int ebN_global_new = elementBoundaryNumbering_global_old2new[ebmp->first];
4788 assert(elementBoundaryNumbering_global2subdomainMap.find(ebN_global_new) != elementBoundaryNumbering_global2subdomainMap.end());
4789 int ebN_subdomain = elementBoundaryNumbering_global2subdomainMap[ebN_global_new];
4792 if (!hasVertexMarkers)
4803 elementBoundaryNumbering_global_old2new.resize(0);
4805 PetscLogEventEnd(build_subdomains_renumber_event,0,0,0,0);
4806 int build_subdomains_cleanup_event;
4807 PetscLogEventRegister(
"Cleanup",0,&build_subdomains_cleanup_event);
4808 PetscLogEventBegin(build_subdomains_cleanup_event,0,0,0,0);
4822 for (
int sdN = 0; sdN < size+1; sdN++)
4853 H5Sclose(filespace);
4858 PetscLogEventEnd(build_subdomains_cleanup_event,0,0,0,0);
4860 PetscLogView(PETSC_VIEWER_STDOUT_WORLD);
4867 using namespace std;
4870 ierr = MPI_Comm_size(PROTEUS_COMM_WORLD,&size);
4871 ierr = MPI_Comm_rank(PROTEUS_COMM_WORLD,&rank);
4905 valarray<int> elementOffsets_old(size+1);
4906 elementOffsets_old[0] = 0;
4907 for(
int sdN=0;sdN<size;sdN++)
4908 elementOffsets_old[sdN+1] = elementOffsets_old[sdN] +
4914 int nElements_subdomain = (elementOffsets_old[rank+1] - elementOffsets_old[rank]);
4915 PetscInt *elementNeighborsOffsets_subdomain,*elementNeighbors_subdomain,*weights_subdomain;
4916 PetscMalloc(
sizeof(PetscInt)*(nElements_subdomain+1),&elementNeighborsOffsets_subdomain);
4920 elementNeighborsOffsets_subdomain[0] = 0;
4921 for (
int eN=0,offset=0; eN < nElements_subdomain; eN++)
4923 int eN_global = elementOffsets_old[rank] + eN;
4924 int offsetStart=offset;
4928 if (eN_neighbor_global >= 0 )
4929 elementNeighbors_subdomain[offset++] = eN_neighbor_global;
4931 elementNeighborsOffsets_subdomain[eN+1]=offset;
4932 sort(&elementNeighbors_subdomain[offsetStart],&elementNeighbors_subdomain[offset]);
4933 int weight = (elementNeighborsOffsets_subdomain[eN+1] - elementNeighborsOffsets_subdomain[eN]);
4934 for (
int k=elementNeighborsOffsets_subdomain[eN];k<elementNeighborsOffsets_subdomain[eN+1];k++)
4935 weights_subdomain[k] = weight;
4944 ierr = MatCreateMPIAdj(PROTEUS_COMM_WORLD,
4945 nElements_subdomain,
4947 elementNeighborsOffsets_subdomain,
4948 elementNeighbors_subdomain,
4950 &petscAdjacency);CHKERRABORT(PROTEUS_COMM_WORLD, ierr);
4951 PetscFree(weights_subdomain);
4952 MatPartitioning petscPartition;
4953 MatPartitioningCreate(PROTEUS_COMM_WORLD,&petscPartition);
4954 MatPartitioningSetAdjacency(petscPartition,petscAdjacency);
4955 MatPartitioningSetFromOptions(petscPartition);
4958 IS elementPartitioningIS_new;
4959 MatPartitioningApply(petscPartition,&elementPartitioningIS_new);
4960 MatPartitioningDestroy(&petscPartition);
4961 MatDestroy(&petscAdjacency);
4993 valarray<int> nElements_subdomain_new(size);
4994 ISPartitioningCount(elementPartitioningIS_new,size,&nElements_subdomain_new[0]);
4997 valarray<int> elementOffsets_new(size+1);
4998 elementOffsets_new[0] = 0;
4999 for (
int sdN=0;sdN<size;sdN++)
5000 elementOffsets_new[sdN+1] = elementOffsets_new[sdN] + nElements_subdomain_new[sdN];
5003 IS elementNumberingIS_subdomain_old2new;
5004 ISPartitioningToNumbering(elementPartitioningIS_new,&elementNumberingIS_subdomain_old2new);
5011 IS elementNumberingIS_global_old2new;
5012 ISAllGather(elementNumberingIS_subdomain_old2new,&elementNumberingIS_global_old2new);
5013 const PetscInt *elementNumbering_global_old2new;
5014 ISGetIndices(elementNumberingIS_global_old2new,&elementNumbering_global_old2new);
5017 elementNumbering_global_new2old[elementNumbering_global_old2new[eN]] = eN;
5055 elementNumbering_global_old2new[eN_ebN];
5067 MPI_Recv(nodeMask,PetscBTLength(mesh.
nNodes_global),MPI_CHAR,rank-1,0,PROTEUS_COMM_WORLD,&status);
5070 set<int> nodes_subdomain_owned;
5071 for(
int eN=elementOffsets_new[rank];eN<elementOffsets_new[rank+1];eN++)
5074 int nN_global = elementNodesArray_new[eN*mesh.
nNodes_element+nN];
5075 if (!PetscBTLookupSet(nodeMask,nN_global))
5076 nodes_subdomain_owned.insert(nN_global);
5080 MPI_Send(nodeMask,PetscBTLength(mesh.
nNodes_global),MPI_CHAR,rank+1,0,PROTEUS_COMM_WORLD);
5081 ierr = PetscBTDestroy(&nodeMask);
5083 cerr<<
"Error in PetscBTDestroy"<<endl;
5085 valarray<int> nNodes_subdomain_new(size),
5086 nodeOffsets_new(size+1);
5087 for (
int sdN=0;sdN<size;sdN++)
5089 nNodes_subdomain_new[sdN] = nodes_subdomain_owned.size();
5091 nNodes_subdomain_new[sdN] = 0;
5092 valarray<int> nNodes_subdomain_new_send=nNodes_subdomain_new;
5093 MPI_Allreduce(&nNodes_subdomain_new_send[0],&nNodes_subdomain_new[0],size,MPI_INT,MPI_SUM,PROTEUS_COMM_WORLD);
5094 nodeOffsets_new[0] = 0;
5095 for (
int sdN=0;sdN<size;sdN++)
5096 nodeOffsets_new[sdN+1] = nodeOffsets_new[sdN]+nNodes_subdomain_new[sdN];
5102 valarray<int> nodeNumbering_new2old(nodes_subdomain_owned.size());
5103 set<int>::iterator nN_ownedp=nodes_subdomain_owned.begin();
5104 for (
int nN=0;nN<int(nodes_subdomain_owned.size());nN++)
5106 nodeNumbering_new2old[nN] = *nN_ownedp++;
5108 IS nodeNumberingIS_new2old;
5109 ISCreateGeneral(PROTEUS_COMM_WORLD,nodes_subdomain_owned.size(),&nodeNumbering_new2old[0],PETSC_COPY_VALUES,&nodeNumberingIS_new2old);
5110 IS nodeNumberingIS_global_new2old;
5111 ISAllGather(nodeNumberingIS_new2old,&nodeNumberingIS_global_new2old);
5112 const PetscInt *nodeNumbering_global_new2old;
5113 valarray<int> nodeNumbering_old2new_global(mesh.
nNodes_global);
5114 ISGetIndices(nodeNumberingIS_global_new2old,&nodeNumbering_global_new2old);
5117 nodeNumbering_old2new_global[nodeNumbering_global_new2old[nN]] = nN;
5125 elementNodesArray_new[eN*mesh.
nNodes_element+nN] = nodeNumbering_old2new_global[nN_old];
5137 edgeNodesArray_new[i] = nodeNumbering_old2new_global[nN_old];
5143 nodeStarArray_new[i] = nodeNumbering_old2new_global[nN_old];
5149 int nN_new = nodeNumbering_old2new_global[nN];
5150 nodeArray_new[nN_new*3+0] = mesh.
nodeArray[nN*3+0];
5151 nodeArray_new[nN_new*3+1] = mesh.
nodeArray[nN*3+1];
5152 nodeArray_new[nN_new*3+2] = mesh.
nodeArray[nN*3+2];
5157 MPI_Status status_elementBoundaries;
5158 PetscBT elementBoundaryMask;
5162 MPI_Recv(elementBoundaryMask,PetscBTLength(mesh.
nElementBoundaries_global),MPI_CHAR,rank-1,0,PROTEUS_COMM_WORLD,&status_elementBoundaries);
5165 set<int> elementBoundaries_subdomain_owned;
5166 for(
int eN=elementOffsets_new[rank];eN<elementOffsets_new[rank+1];eN++)
5170 if(!PetscBTLookup(elementBoundaryMask,ebN_global))
5176 if (nodes_subdomain_owned.count(nN_global_old) > 0)
5184 PetscBTSet(elementBoundaryMask,ebN_global);
5185 elementBoundaries_subdomain_owned.insert(ebN_global);
5189 PetscBTSet(elementBoundaryMask,ebN_global);
5190 elementBoundaries_subdomain_owned.insert(ebN_global);
5198 ierr = PetscBTDestroy(&elementBoundaryMask);
5200 cerr<<
"Error in PetscBTDestroy for elementBoundaries"<<endl;
5202 valarray<int> nElementBoundaries_subdomain_new(size),
5203 elementBoundaryOffsets_new(size+1);
5204 for (
int sdN=0;sdN<size;sdN++)
5206 nElementBoundaries_subdomain_new[sdN] = elementBoundaries_subdomain_owned.size();
5208 nElementBoundaries_subdomain_new[sdN] = 0;
5209 valarray<int> nElementBoundaries_subdomain_new_send=nElementBoundaries_subdomain_new;
5210 MPI_Allreduce(&nElementBoundaries_subdomain_new_send[0],&nElementBoundaries_subdomain_new[0],size,MPI_INT,MPI_SUM,PROTEUS_COMM_WORLD);
5211 elementBoundaryOffsets_new[0] = 0;
5212 for (
int sdN=0;sdN<size;sdN++)
5213 elementBoundaryOffsets_new[sdN+1] = elementBoundaryOffsets_new[sdN]+nElementBoundaries_subdomain_new[sdN];
5218 valarray<int> elementBoundaryNumbering_new2old(elementBoundaries_subdomain_owned.size());
5219 set<int>::iterator ebN_ownedp=elementBoundaries_subdomain_owned.begin();
5220 for (
int ebN=0;ebN<int(elementBoundaries_subdomain_owned.size());ebN++)
5222 elementBoundaryNumbering_new2old[ebN] = *ebN_ownedp++;
5224 IS elementBoundaryNumberingIS_new2old;
5225 ISCreateGeneral(PROTEUS_COMM_WORLD,elementBoundaries_subdomain_owned.size(),&elementBoundaryNumbering_new2old[0],PETSC_COPY_VALUES,&elementBoundaryNumberingIS_new2old);
5226 IS elementBoundaryNumberingIS_global_new2old;
5227 ISAllGather(elementBoundaryNumberingIS_new2old,&elementBoundaryNumberingIS_global_new2old);
5228 const PetscInt *elementBoundaryNumbering_global_new2old;
5230 ISGetIndices(elementBoundaryNumberingIS_global_new2old,&elementBoundaryNumbering_global_new2old);
5233 elementBoundaryNumbering_old2new_global[elementBoundaryNumbering_global_new2old[ebN]] = ebN;
5247 int ebN_old=elementBoundaryNumbering_global_new2old[ebN];
5255 elementBoundaryElementsArray_new[ebN*2+0] = elementNumbering_global_old2new[eN_L_old];
5257 elementBoundaryElementsArray_new[ebN*2+1] = elementNumbering_global_old2new[eN_R_old];
5318 MPI_Status status_edges;
5323 MPI_Recv(edgesMask,PetscBTLength(mesh.
nEdges_global),MPI_CHAR,rank-1,0,PROTEUS_COMM_WORLD,&status_edges);
5326 map<NodeTuple<2>,
int> nodesEdgeMap_global;
5327 set<int> edges_subdomain_owned;
5331 nodes[0] = edgeNodesArray_new[2*ig];
5332 nodes[1] = edgeNodesArray_new[2*ig+1];
5334 nodesEdgeMap_global[et] = ig;
5336 for(
int eN=elementOffsets_new[rank];eN<elementOffsets_new[rank+1];eN++)
5345 if (nodesEdgeMap_global.find(et) != nodesEdgeMap_global.end())
5347 int edge_global = nodesEdgeMap_global[et];
5348 if(!PetscBTLookup(edgesMask,edge_global))
5350 PetscBTSet(edgesMask,edge_global);
5351 edges_subdomain_owned.insert(nodesEdgeMap_global[et]);
5357 MPI_Send(edgesMask,PetscBTLength(mesh.
nEdges_global),MPI_CHAR,rank+1,0,PROTEUS_COMM_WORLD);
5358 ierr = PetscBTDestroy(&edgesMask);
5360 cerr<<
"Error in PetscBTDestroy for edges"<<endl;
5362 valarray<int> nEdges_subdomain_new(size),
5363 edgeOffsets_new(size+1);
5365 for (
int sdN=0; sdN < size; sdN++)
5368 nEdges_subdomain_new[sdN] = edges_subdomain_owned.size();
5370 nEdges_subdomain_new[sdN] = 0;
5373 valarray<int> nEdges_subdomain_new_send=nEdges_subdomain_new;
5374 MPI_Allreduce(&nEdges_subdomain_new_send[0],&nEdges_subdomain_new[0],size,MPI_INT,MPI_SUM,PROTEUS_COMM_WORLD);
5375 edgeOffsets_new[0] = 0;
5376 for (
int sdN=0;sdN<size;sdN++)
5377 edgeOffsets_new[sdN+1] = edgeOffsets_new[sdN]+nEdges_subdomain_new[sdN];
5380 valarray<int> edgeNumbering_new2old(edges_subdomain_owned.size());
5381 set<int>::iterator edges_ownedp = edges_subdomain_owned.begin();
5382 for (
int i=0; i < int(edges_subdomain_owned.size());i++)
5383 edgeNumbering_new2old[i] = *edges_ownedp++;
5385 IS edgeNumberingIS_new2old;
5386 ISCreateGeneral(PROTEUS_COMM_WORLD,edges_subdomain_owned.size(),&edgeNumbering_new2old[0],PETSC_COPY_VALUES,&edgeNumberingIS_new2old);
5387 IS edgeNumberingIS_global_new2old;
5388 ISAllGather(edgeNumberingIS_new2old,&edgeNumberingIS_global_new2old);
5389 const PetscInt *edgeNumbering_global_new2old;
5392 valarray<int> edgeNumbering_old2new_global(mesh.
nEdges_global);
5393 ISGetIndices(edgeNumberingIS_global_new2old,&edgeNumbering_global_new2old);
5396 edgeNumbering_old2new_global[edgeNumbering_global_new2old[ig]] = ig;
5401 valarray<int> edgeNodesArray_newNodesAndEdges(2*mesh.
nEdges_global);
5402 map<NodeTuple<2>,
int > nodesEdgeMap_global_new;
5406 const int edge_old = edgeNumbering_global_new2old[ig];
5407 nodes[0] = edgeNodesArray_new[edge_old*2+0];
5408 nodes[1] = edgeNodesArray_new[edge_old*2+1];
5410 edgeNodesArray_newNodesAndEdges[ig*2+0] = edgeNodesArray_new[edge_old*2+0];
5411 edgeNodesArray_newNodesAndEdges[ig*2+1] = edgeNodesArray_new[edge_old*2+1];
5412 assert(nodesEdgeMap_global_new.find(et) == nodesEdgeMap_global_new.end());
5413 nodesEdgeMap_global_new[et] = ig;
5418 set<int> elements_overlap,nodes_overlap,elementBoundaries_overlap,edges_overlap;
5419 for(
int eN=elementOffsets_new[rank];eN<elementOffsets_new[rank+1];eN++)
5423 int nN_global = elementNodesArray_new[eN*mesh.
nNodes_element+nN];
5424 if (nN_global < nodeOffsets_new[rank] || nN_global >= nodeOffsets_new[rank+1])
5425 nodes_overlap.insert(nN_global);
5430 if (ebN_global < elementBoundaryOffsets_new[rank] || ebN_global >= elementBoundaryOffsets_new[rank+1])
5431 elementBoundaries_overlap.insert(ebN_global);
5440 const int nN0_global = elementNodesArray_new[eN*mesh.
nNodes_element+nN0];
5441 const int nN1_global = elementNodesArray_new[eN*mesh.
nNodes_element+nN1];
5442 if (nN0_global < nodeOffsets_new[rank] || nN0_global >= nodeOffsets_new[rank+1])
5443 assert(nodes_overlap.find(nN0_global) != nodes_overlap.end());
5444 if (nN1_global < nodeOffsets_new[rank] || nN1_global >= nodeOffsets_new[rank+1])
5445 assert(nodes_overlap.find(nN1_global) != nodes_overlap.end());
5446 if (nodesEdgeMap_global_new.find(et) != nodesEdgeMap_global_new.end())
5448 const int edge_global = nodesEdgeMap_global_new[et];
5449 if (edge_global < edgeOffsets_new[rank] || edge_global >= edgeOffsets_new[rank+1])
5450 edges_overlap.insert(edge_global);
5455 if (nElements_overlap > 0)
5478 for (set<int>::iterator nN=nodes_subdomain_owned.begin();nN != nodes_subdomain_owned.end();nN++)
5483 int eN_new = elementNumbering_global_old2new[eN];
5484 if (eN_new < elementOffsets_new[rank] or eN_new >= elementOffsets_new[rank+1])
5486 elements_overlap.insert(eN_new);
5487 for (
int nN_element=0;nN_element<mesh.
nNodes_element;nN_element++)
5489 int nN_global = elementNodesArray_new[eN_new*mesh.
nNodes_element+nN_element];
5490 if (nN_global < nodeOffsets_new[rank] or nN_global >= nodeOffsets_new[rank+1])
5491 nodes_overlap.insert(nN_global);
5496 if (ebN_global < elementBoundaryOffsets_new[rank] || ebN_global >= elementBoundaryOffsets_new[rank+1])
5497 elementBoundaries_overlap.insert(ebN_global);
5504 nodes[0] = elementNodesArray_new[eN_new*mesh.
nNodes_element+nN0];
5505 nodes[1] = elementNodesArray_new[eN_new*mesh.
nNodes_element+nN1];
5507 const int nN0_global = elementNodesArray_new[eN_new*mesh.
nNodes_element+nN0];
5508 const int nN1_global = elementNodesArray_new[eN_new*mesh.
nNodes_element+nN1];
5509 if (nodesEdgeMap_global_new.find(et) != nodesEdgeMap_global_new.end())
5511 const int edge_global = nodesEdgeMap_global_new[et];
5512 if (edge_global < edgeOffsets_new[rank] || edge_global >= edgeOffsets_new[rank+1])
5513 edges_overlap.insert(edge_global);
5520 for(
int eN=elementOffsets_new[rank];eN<elementOffsets_new[rank+1];eN++)
5525 (eN_ebN < elementOffsets_new[rank] || eN_ebN >= elementOffsets_new[rank+1]))
5527 elements_overlap.insert(eN_ebN);
5530 int nN_global = elementNodesArray_new[eN_ebN*mesh.
nNodes_element+nN];
5531 if (nN_global < nodeOffsets_new[rank] || nN_global >= nodeOffsets_new[rank+1])
5532 nodes_overlap.insert(nN_global);
5537 if (ebN_global < elementBoundaryOffsets_new[rank] || ebN_global >= elementBoundaryOffsets_new[rank+1])
5538 elementBoundaries_overlap.insert(ebN_global);
5545 nodes[0] = elementNodesArray_new[eN_ebN*mesh.
nNodes_element+nN0];
5546 nodes[1] = elementNodesArray_new[eN_ebN*mesh.
nNodes_element+nN1];
5548 const int nN0_global = elementNodesArray_new[eN_ebN*mesh.
nNodes_element+nN0];
5549 const int nN1_global = elementNodesArray_new[eN_ebN*mesh.
nNodes_element+nN1];
5550 bool foundEdge =
false;
5551 const int edge_global = nodesEdgeMap_global_new[et];
5552 if (edge_global < edgeOffsets_new[rank] || edge_global >= edgeOffsets_new[rank+1])
5553 edges_overlap.insert(edge_global);
5558 for (
int layer=1;layer<nElements_overlap;layer++)
5560 for (set<int>::iterator eN_p=elements_overlap.begin();eN_p != elements_overlap.end();eN_p++)
5562 int eN_global = *eN_p;
5567 (eN_ebN < elementOffsets_new[rank] || eN_ebN >= elementOffsets_new[rank+1]))
5569 elements_overlap.insert(eN_ebN);
5572 int nN_global = elementNodesArray_new[eN_ebN*mesh.
nNodes_element+nN];
5573 if (nN_global < nodeOffsets_new[rank] || nN_global >= nodeOffsets_new[rank+1])
5574 nodes_overlap.insert(nN_global);
5579 if (ebN_global < elementBoundaryOffsets_new[rank] || ebN_global >= elementBoundaryOffsets_new[rank+1])
5580 elementBoundaries_overlap.insert(ebN_global);
5587 nodes[0] = elementNodesArray_new[eN_ebN*mesh.
nNodes_element+nN0];
5588 nodes[1] = elementNodesArray_new[eN_ebN*mesh.
nNodes_element+nN1];
5590 const int nN0_global = elementNodesArray_new[eN_ebN*mesh.
nNodes_element+nN0];
5591 const int nN1_global = elementNodesArray_new[eN_ebN*mesh.
nNodes_element+nN1];
5592 bool foundEdge =
false;
5593 if (nodesEdgeMap_global_new.find(et) != nodesEdgeMap_global_new.end())
5595 const int edge_global = nodesEdgeMap_global_new[et];
5596 if (edge_global < edgeOffsets_new[rank] || edge_global >= edgeOffsets_new[rank+1])
5597 edges_overlap.insert(edge_global);
5621 map<int,int> nodeNumbering_global2subdomain;
5624 for(
int nN=0;nN<nNodes_subdomain_new[rank];nN++)
5626 int nN_global = nN + nodeOffsets_new[rank];
5627 nodeNumbering_subdomain2global[nN] = nN_global;
5628 nodeNumbering_global2subdomain[nN_global] = nN;
5636 set<int>::iterator nN_p=nodes_overlap.begin();
5637 for(
int nN=nNodes_subdomain_new[rank];nN < nNodes_subdomain_new[rank] + int(nodes_overlap.size()); nN++)
5639 int nN_global = *nN_p++;
5640 nodeNumbering_subdomain2global[nN] = nN_global;
5641 nodeNumbering_global2subdomain[nN_global] = nN;
5650 map<int,int> elementBoundaryNumbering_global2subdomain;
5651 for (
int ebN=0; ebN < nElementBoundaries_subdomain_new[rank]; ebN++)
5653 int ebN_global = ebN + elementBoundaryOffsets_new[rank];
5654 elementBoundaryNumbering_subdomain2global[ebN]=ebN_global;
5655 elementBoundaryNumbering_global2subdomain[ebN_global] = ebN;
5657 set<int>::iterator ebN_p = elementBoundaries_overlap.begin();
5658 for(
int ebN=nElementBoundaries_subdomain_new[rank];ebN < nElementBoundaries_subdomain_new[rank] + int(elementBoundaries_overlap.size()); ebN++)
5660 int ebN_global = *ebN_p++;
5661 elementBoundaryNumbering_subdomain2global[ebN] = ebN_global;
5662 elementBoundaryNumbering_global2subdomain[ebN_global] = ebN;
5672 for (
int eN=0;eN<nElements_subdomain_new[rank];eN++)
5674 int eN_global = eN+elementOffsets_new[rank];
5675 elementNumbering_subdomain2global[eN] = eN_global;
5680 nodeNumbering_global2subdomain[elementNodesArray_new[eN_global*mesh.
nNodes_element + nN]];
5686 set<int>::iterator eN_p=elements_overlap.begin();
5687 for(
int eN=nElements_subdomain_new[rank];eN < nElements_subdomain_new[rank]+int(elements_overlap.size());eN++)
5689 int eN_global = *eN_p++;
5691 elementNumbering_subdomain2global[eN] = eN_global;
5695 nodeNumbering_global2subdomain[elementNodesArray_new[eN_global*mesh.
nNodes_element + nN]];
5706 for (
int i=0; i < nEdges_subdomain_new[rank]; i++)
5708 const int ig = i+edgeOffsets_new[rank];
5709 const int nN0_global = edgeNodesArray_newNodesAndEdges[ig*2+0];
5710 const int nN1_global = edgeNodesArray_newNodesAndEdges[ig*2+1];
5712 assert(nodeNumbering_global2subdomain.find(nN0_global) != nodeNumbering_global2subdomain.end());
5713 assert(nodeNumbering_global2subdomain.find(nN1_global) != nodeNumbering_global2subdomain.end());
5714 const int nN0_subdomain = nodeNumbering_global2subdomain[nN0_global];
5715 const int nN1_subdomain = nodeNumbering_global2subdomain[nN1_global];
5718 edgeNumbering_subdomain2global[i] = ig;
5720 set<int>::iterator edge_p = edges_overlap.begin();
5721 for (
int i=nEdges_subdomain_new[rank]; i < nEdges_subdomain_new[rank] + int(edges_overlap.size()); i++)
5723 const int ig =*edge_p++;
5724 const int nN0_global = edgeNodesArray_newNodesAndEdges[ig*2+0];
5725 const int nN1_global = edgeNodesArray_newNodesAndEdges[ig*2+1];
5727 const int nN0_subdomain = nodeNumbering_global2subdomain[nN0_global];
5728 const int nN1_subdomain = nodeNumbering_global2subdomain[nN1_global];
5731 edgeNumbering_subdomain2global[i] = ig;
5792 int eN_global_new = elementNumbering_subdomain2global[eN];
5793 int eN_global_old = elementNumbering_global_new2old[eN_global_new];
5809 for (
int sdN=0;sdN<size+1;sdN++)
5876 ISRestoreIndices(elementNumberingIS_global_old2new,&elementNumbering_global_old2new);
5878 ISDestroy(&elementPartitioningIS_new);
5879 ISDestroy(&elementNumberingIS_subdomain_old2new);
5880 ISDestroy(&elementNumberingIS_global_old2new);
5882 ISRestoreIndices(nodeNumberingIS_global_new2old,&nodeNumbering_global_new2old);
5884 ISDestroy(&nodeNumberingIS_new2old);
5885 ISDestroy(&nodeNumberingIS_global_new2old);
5887 ISRestoreIndices(elementBoundaryNumberingIS_global_new2old,&elementBoundaryNumbering_global_new2old);
5889 ISDestroy(&elementBoundaryNumberingIS_new2old);
5890 ISDestroy(&elementBoundaryNumberingIS_global_new2old);
5892 ISRestoreIndices(edgeNumberingIS_global_new2old,&edgeNumbering_global_new2old);
5894 ISDestroy(&edgeNumberingIS_new2old);
5895 ISDestroy(&edgeNumberingIS_global_new2old);