11 ierr = MPI_Comm_size(PROTEUS_COMM_WORLD,&size);
12 ierr = MPI_Comm_rank(PROTEUS_COMM_WORLD,&rank);
46 valarray<int> elementOffsets_old(size+1);
47 elementOffsets_old[0] = 0;
48 for(
int sdN=0;sdN<size;sdN++)
49 elementOffsets_old[sdN+1] = elementOffsets_old[sdN] +
55 int nElements_subdomain = (elementOffsets_old[rank+1] - elementOffsets_old[rank]);
56 PetscInt *elementNeighborsOffsets_subdomain,*elementNeighbors_subdomain,*weights_subdomain;
57 PetscMalloc(
sizeof(PetscInt)*(nElements_subdomain+1),&elementNeighborsOffsets_subdomain);
61 elementNeighborsOffsets_subdomain[0] = 0;
62 for (
int eN=0,offset=0; eN < nElements_subdomain; eN++)
64 int eN_global = elementOffsets_old[rank] + eN;
65 int offsetStart=offset;
69 if (eN_neighbor_global >= 0 )
70 elementNeighbors_subdomain[offset++] = eN_neighbor_global;
72 elementNeighborsOffsets_subdomain[eN+1]=offset;
73 sort(&elementNeighbors_subdomain[offsetStart],&elementNeighbors_subdomain[offset]);
74 int weight = (elementNeighborsOffsets_subdomain[eN+1] - elementNeighborsOffsets_subdomain[eN]);
75 for (
int k=elementNeighborsOffsets_subdomain[eN];k < elementNeighborsOffsets_subdomain[eN+1];k++)
76 weights_subdomain[k] = weight;
90 ierr = MatCreateMPIAdj(PROTEUS_COMM_WORLD,
93 elementNeighborsOffsets_subdomain,
94 elementNeighbors_subdomain,
96 &petscAdjacency);CHKERRABORT(PROTEUS_COMM_WORLD, ierr);
97 MatPartitioning petscPartition;
98 MatPartitioningCreate(PROTEUS_COMM_WORLD,&petscPartition);
99 MatPartitioningSetAdjacency(petscPartition,petscAdjacency);
100 MatPartitioningSetFromOptions(petscPartition);
103 IS elementPartitioningIS_new;
104 MatPartitioningApply(petscPartition,&elementPartitioningIS_new);
105 MatPartitioningDestroy(&petscPartition);
138 valarray<int> nElements_subdomain_new(size);
139 ISPartitioningCount(elementPartitioningIS_new,size,&nElements_subdomain_new[0]);
142 valarray<int> elementOffsets_new(size+1);
143 elementOffsets_new[0] = 0;
144 for (
int sdN=0;sdN<size;sdN++)
145 elementOffsets_new[sdN+1] = elementOffsets_new[sdN] + nElements_subdomain_new[sdN];
148 IS elementNumberingIS_subdomain_old2new;
149 ISPartitioningToNumbering(elementPartitioningIS_new,&elementNumberingIS_subdomain_old2new);
156 IS elementNumberingIS_global_old2new;
157 ISAllGather(elementNumberingIS_subdomain_old2new,&elementNumberingIS_global_old2new);
159 const PetscInt *elementNumbering_global_old2new;
160 ISGetIndices(elementNumberingIS_global_old2new,&elementNumbering_global_old2new);
163 elementNumbering_global_new2old[elementNumbering_global_old2new[eN]] = eN;
193 elementNumbering_global_old2new[eN_ebN];
200 elementBoundaryElementsArray_new[ebN*2+0] = elementNumbering_global_old2new[eN_L_old];
202 elementBoundaryElementsArray_new[ebN*2+1] = elementNumbering_global_old2new[eN_R_old];
215 MPI_Recv(nodeMask,PetscBTLength(mesh.
nNodes_global),MPI_CHAR,rank-1,0,PROTEUS_COMM_WORLD,&status);
218 set<int> nodes_subdomain_owned;
219 for(
int eN=elementOffsets_new[rank];eN<elementOffsets_new[rank+1];eN++)
223 if (!PetscBTLookupSet(nodeMask,nN_global))
224 nodes_subdomain_owned.insert(nN_global);
228 MPI_Send(nodeMask,PetscBTLength(mesh.
nNodes_global),MPI_CHAR,rank+1,0,PROTEUS_COMM_WORLD);
229 ierr = PetscBTDestroy(&nodeMask);
231 cerr<<
"Error in PetscBTDestroy"<<endl;
233 valarray<int> nNodes_subdomain_new(size),
234 nodeOffsets_new(size+1);
235 for (
int sdN=0;sdN<size;sdN++)
237 nNodes_subdomain_new[sdN] = nodes_subdomain_owned.size();
239 nNodes_subdomain_new[sdN] = 0;
240 valarray<int> nNodes_subdomain_new_send=nNodes_subdomain_new;
241 MPI_Allreduce(&nNodes_subdomain_new_send[0],&nNodes_subdomain_new[0],size,MPI_INT,MPI_SUM,PROTEUS_COMM_WORLD);
242 nodeOffsets_new[0] = 0;
243 for (
int sdN=0;sdN<size;sdN++)
244 nodeOffsets_new[sdN+1] = nodeOffsets_new[sdN]+nNodes_subdomain_new[sdN];
247 valarray<int> nodeNumbering_new2old(nodes_subdomain_owned.size());
248 set<int>::iterator nN_ownedp=nodes_subdomain_owned.begin();
249 for (
int nN=0;nN<int(nodes_subdomain_owned.size());nN++)
251 nodeNumbering_new2old[nN] = *nN_ownedp++;
253 IS nodeNumberingIS_new2old;
254 ISCreateGeneral(PROTEUS_COMM_WORLD,nodes_subdomain_owned.size(),&nodeNumbering_new2old[0],PETSC_COPY_VALUES,&nodeNumberingIS_new2old);
255 IS nodeNumberingIS_global_new2old;
256 ISAllGather(nodeNumberingIS_new2old,&nodeNumberingIS_global_new2old);
257 const PetscInt *nodeNumbering_global_new2old;
258 valarray<int> nodeNumbering_old2new_global(mesh.
nNodes_global);
259 ISGetIndices(nodeNumberingIS_global_new2old,&nodeNumbering_global_new2old);
262 nodeNumbering_old2new_global[nodeNumbering_global_new2old[nN]] = nN;
270 elementNodesArray_new[eN*mesh.
nNodes_element+nN] = nodeNumbering_old2new_global[nN_old];
276 elementBoundaryNodesArray_new[i] = nodeNumbering_old2new_global[nN_old];
281 edgeNodesArray_new[i] = nodeNumbering_old2new_global[nN_old];
287 nodeStarArray_new[i] = nodeNumbering_old2new_global[nN_old];
293 int nN_new = nodeNumbering_old2new_global[nN];
294 nodeArray_new[nN_new*3+0] = mesh.
nodeArray[nN*3+0];
295 nodeArray_new[nN_new*3+1] = mesh.
nodeArray[nN*3+1];
296 nodeArray_new[nN_new*3+2] = mesh.
nodeArray[nN*3+2];
321 set<int> elements_overlap,nodes_overlap;
322 for(
int eN=elementOffsets_new[rank];eN<elementOffsets_new[rank+1];eN++)
326 if (nN_global < nodeOffsets_new[rank] || nN_global >= nodeOffsets_new[rank+1])
327 nodes_overlap.insert(nN_global);
331 if (nElements_overlap > 0)
354 for (set<int>::iterator nN=nodes_subdomain_owned.begin();nN != nodes_subdomain_owned.end();nN++)
359 int eN_new = elementNumbering_global_old2new[eN];
360 if (eN_new < elementOffsets_new[rank] or eN_new >= elementOffsets_new[rank+1])
362 elements_overlap.insert(eN_new);
363 for (
int nN_element=0;nN_element<mesh.
nNodes_element;nN_element++)
365 int nN_global = elementNodesArray_new[eN_new*mesh.
nNodes_element+nN_element];
366 if (nN_global < nodeOffsets_new[rank] or nN_global >= nodeOffsets_new[rank+1])
367 nodes_overlap.insert(nN_global);
373 for(
int eN=elementOffsets_new[rank];eN<elementOffsets_new[rank+1];eN++)
378 (eN_ebN < elementOffsets_new[rank] || eN_ebN >= elementOffsets_new[rank+1]))
380 elements_overlap.insert(eN_ebN);
383 int nN_global = elementNodesArray_new[eN_ebN*mesh.
nNodes_element+nN];
384 if (nN_global < nodeOffsets_new[rank] || nN_global >= nodeOffsets_new[rank+1])
385 nodes_overlap.insert(nN_global);
390 for (
int layer=1;layer<nElements_overlap;layer++)
392 for (set<int>::iterator eN_p=elements_overlap.begin();eN_p != elements_overlap.end();eN_p++)
394 int eN_global = *eN_p;
399 (eN_ebN < elementOffsets_new[rank] || eN_ebN >= elementOffsets_new[rank+1]))
401 elements_overlap.insert(eN_ebN);
404 int nN_global = elementNodesArray_new[eN_ebN*mesh.
nNodes_element+nN];
405 if (nN_global < nodeOffsets_new[rank] || nN_global >= nodeOffsets_new[rank+1])
406 nodes_overlap.insert(nN_global);
424 map<int,int> nodeNumbering_global2subdomain;
427 for(
int nN=0;nN<nNodes_subdomain_new[rank];nN++)
429 int nN_global = nN + nodeOffsets_new[rank];
430 nodeNumbering_subdomain2global[nN] = nN_global;
431 nodeNumbering_global2subdomain[nN_global] = nN;
439 set<int>::iterator nN_p=nodes_overlap.begin();
440 for(
int nN=nNodes_subdomain_new[rank];nN < nNodes_subdomain_new[rank] + int(nodes_overlap.size()); nN++)
442 int nN_global = *nN_p++;
443 nodeNumbering_subdomain2global[nN] = nN_global;
444 nodeNumbering_global2subdomain[nN_global] = nN;
453 for (
int eN=0;eN<nElements_subdomain_new[rank];eN++)
455 int eN_global = eN+elementOffsets_new[rank];
456 elementNumbering_subdomain2global[eN] = eN_global;
460 nodeNumbering_global2subdomain[elementNodesArray_new[eN_global*mesh.
nNodes_element + nN]];
462 set<int>::iterator eN_p=elements_overlap.begin();
463 for(
int eN=nElements_subdomain_new[rank];eN < nElements_subdomain_new[rank]+int(elements_overlap.size());eN++)
465 int eN_global = *eN_p++;
467 elementNumbering_subdomain2global[eN] = eN_global;
470 nodeNumbering_global2subdomain[elementNodesArray_new[eN_global*mesh.
nNodes_element + nN]];
497 int eN_global_new = elementNumbering_subdomain2global[eN];
498 int eN_global_old = elementNumbering_global_new2old[eN_global_new];
512 for (
int sdN=0;sdN<size+1;sdN++)
525 ISRestoreIndices(elementNumberingIS_global_old2new,&elementNumbering_global_old2new);
527 ISDestroy(&elementPartitioningIS_new);
528 ISDestroy(&elementNumberingIS_subdomain_old2new);
529 ISDestroy(&elementNumberingIS_global_old2new);
531 ISRestoreIndices(nodeNumberingIS_global_new2old,&nodeNumbering_global_new2old);
533 ISDestroy(&nodeNumberingIS_new2old);
534 ISDestroy(&nodeNumberingIS_global_new2old);
544 ierr = MPI_Comm_size(PROTEUS_COMM_WORLD,&size);
545 ierr = MPI_Comm_rank(PROTEUS_COMM_WORLD,&rank);
574 valarray<int> nodeOffsets_old(size+1);
575 nodeOffsets_old[0] = 0;
576 for (
int sdN=0; sdN < size; sdN++)
578 nodeOffsets_old[sdN+1] = nodeOffsets_old[sdN] +
584 int nNodes_subdomain = (nodeOffsets_old[rank+1] - nodeOffsets_old[rank]);
585 PetscInt *nodeNeighborsOffsets_subdomain,*nodeNeighbors_subdomain,*weights_subdomain;
586 PetscMalloc(
sizeof(PetscInt)*(nNodes_subdomain+1),&nodeNeighborsOffsets_subdomain);
589 nodeNeighborsOffsets_subdomain[0] = 0;
590 for (
int nN = 0,offset=0; nN < nNodes_subdomain; nN++)
592 int nN_global = nodeOffsets_old[rank] + nN;
596 nodeNeighbors_subdomain[offset++] = mesh.
nodeStarArray[offset_global];
598 nodeNeighborsOffsets_subdomain[nN+1]=offset;
599 sort(&nodeNeighbors_subdomain[nodeNeighborsOffsets_subdomain[nN]],&nodeNeighbors_subdomain[nodeNeighborsOffsets_subdomain[nN+1]]);
600 int weight= (nodeNeighborsOffsets_subdomain[nN+1] - nodeNeighborsOffsets_subdomain[nN]);
601 for (
int k=nodeNeighborsOffsets_subdomain[nN];k<nodeNeighborsOffsets_subdomain[nN+1];k++)
602 weights_subdomain[k] = weight;
613 ierr = MatCreateMPIAdj(PROTEUS_COMM_WORLD,
616 nodeNeighborsOffsets_subdomain,
617 nodeNeighbors_subdomain,
619 &petscAdjacency);CHKERRABORT(PROTEUS_COMM_WORLD, ierr);
620 MatPartitioning petscPartition;
621 MatPartitioningCreate(PROTEUS_COMM_WORLD,&petscPartition);
622 MatPartitioningSetAdjacency(petscPartition,petscAdjacency);
623 MatPartitioningSetFromOptions(petscPartition);
626 IS nodePartitioningIS_new;
627 MatPartitioningApply(petscPartition,&nodePartitioningIS_new);
628 MatPartitioningDestroy(&petscPartition);
631 valarray<int> nNodes_subdomain_new(size);
632 ISPartitioningCount(nodePartitioningIS_new,size,&nNodes_subdomain_new[0]);
635 valarray<int> nodeOffsets_new(size+1);
636 nodeOffsets_new[0] = 0;
637 for (
int sdN = 0; sdN < size; sdN++)
638 nodeOffsets_new[sdN+1] = nodeOffsets_new[sdN] + nNodes_subdomain_new[sdN];
641 IS nodeNumberingIS_subdomain_old2new;
642 ISPartitioningToNumbering(nodePartitioningIS_new,&nodeNumberingIS_subdomain_old2new);
647 IS nodeNumberingIS_global_old2new;
648 ISAllGather(nodeNumberingIS_subdomain_old2new,&nodeNumberingIS_global_old2new);
650 const PetscInt * nodeNumbering_global_old2new;
651 ISGetIndices(nodeNumberingIS_global_old2new,&nodeNumbering_global_old2new);
654 valarray<int> nodeNumbering_global_new2old(mesh.
nNodes_global);
656 nodeNumbering_global_new2old[nodeNumbering_global_old2new[nN]] = nN;
691 MPI_Recv(elementMask,PetscBTLength(mesh.
nElements_global),MPI_CHAR,rank-1,0,PROTEUS_COMM_WORLD,&status);
694 set<int> elements_subdomain_owned;
695 for (
int nN = nodeOffsets_new[rank]; nN < nodeOffsets_new[rank+1]; nN++)
697 int nN_global_old = nodeNumbering_global_new2old[nN];
702 if (!PetscBTLookupSet(elementMask,eN_star_old))
704 elements_subdomain_owned.insert(eN_star_old);
735 MPI_Send(elementMask,PetscBTLength(mesh.
nElements_global),MPI_CHAR,rank+1,0,PROTEUS_COMM_WORLD);
736 ierr = PetscBTDestroy(&elementMask);
738 cerr<<
"Error in PetscBTDestroy"<<endl;
743 valarray<int> nElements_subdomain_new(size),
744 elementOffsets_new(size+1);
745 for (
int sdN = 0; sdN < size; sdN++)
748 nElements_subdomain_new[sdN] = int(elements_subdomain_owned.size());
750 nElements_subdomain_new[sdN] = 0;
752 valarray<int> nElements_subdomain_new_send = nElements_subdomain_new;
753 MPI_Allreduce(&nElements_subdomain_new_send[0],&nElements_subdomain_new[0],size,MPI_INT,MPI_SUM,PROTEUS_COMM_WORLD);
755 elementOffsets_new[0] = 0;
756 for (
int sdN = 0; sdN < size; sdN++)
757 elementOffsets_new[sdN+1] = elementOffsets_new[sdN] + nElements_subdomain_new[sdN];
760 valarray<int> elementNumbering_subdomain_new2old(elements_subdomain_owned.size());
761 set<int>::iterator eN_ownedp = elements_subdomain_owned.begin();
762 for (
int eN = 0; eN < int(elements_subdomain_owned.size()); eN++)
764 elementNumbering_subdomain_new2old[eN] = *eN_ownedp++;
767 IS elementNumberingIS_subdomain_new2old;
768 ISCreateGeneral(PROTEUS_COMM_WORLD,elements_subdomain_owned.size(),&elementNumbering_subdomain_new2old[0],PETSC_COPY_VALUES,
769 &elementNumberingIS_subdomain_new2old);
770 IS elementNumberingIS_global_new2old;
771 ISAllGather(elementNumberingIS_subdomain_new2old,&elementNumberingIS_global_new2old);
773 const PetscInt *elementNumbering_global_new2old;
774 ISGetIndices(elementNumberingIS_global_new2old,&elementNumbering_global_new2old);
779 elementNumbering_global_old2new[elementNumbering_global_new2old[eN]] = eN;
793 MPI_Status status_elementBoundaries;
794 PetscBT elementBoundaryMask;
798 MPI_Recv(elementBoundaryMask,PetscBTLength(mesh.
nElementBoundaries_global),MPI_CHAR,rank-1,0,PROTEUS_COMM_WORLD,&status_elementBoundaries);
802 set<int> elementBoundaries_subdomain_owned;
806 int lface[6][4] = {{0,1,2,3},
814 for (
int nN = nodeOffsets_new[rank]; nN < nodeOffsets_new[rank+1]; nN++)
816 int nN_global_old = nodeNumbering_global_new2old[nN];
822 int eN_star_new = elementNumbering_global_old2new[eN_star_old];
828 bool foundNode =
false;
832 if (nN_global_old_across == nN_global_old) foundNode =
true;
837 if (!PetscBTLookupSet(elementBoundaryMask,ebN_global))
838 elementBoundaries_subdomain_owned.insert(ebN_global);
849 for (
int nN = nodeOffsets_new[rank]; nN < nodeOffsets_new[rank+1]; nN++)
851 int nN_global_old = nodeNumbering_global_new2old[nN];
857 int eN_star_new = elementNumbering_global_old2new[eN_star_old];
864 if (nN_global_old_across != nN_global_old)
867 if (!PetscBTLookupSet(elementBoundaryMask,ebN_global))
868 elementBoundaries_subdomain_owned.insert(ebN_global);
878 ierr = PetscBTDestroy(&elementBoundaryMask);
880 cerr<<
"Error in PetscBTDestroy for elementBoundaries"<<endl;
882 valarray<int> nElementBoundaries_subdomain_new(size),
883 elementBoundaryOffsets_new(size+1);
884 for (
int sdN=0;sdN<size;sdN++)
886 nElementBoundaries_subdomain_new[sdN] = elementBoundaries_subdomain_owned.size();
888 nElementBoundaries_subdomain_new[sdN] = 0;
889 valarray<int> nElementBoundaries_subdomain_new_send=nElementBoundaries_subdomain_new;
890 MPI_Allreduce(&nElementBoundaries_subdomain_new_send[0],&nElementBoundaries_subdomain_new[0],size,MPI_INT,MPI_SUM,PROTEUS_COMM_WORLD);
891 elementBoundaryOffsets_new[0] = 0;
892 for (
int sdN=0;sdN<size;sdN++)
893 elementBoundaryOffsets_new[sdN+1] = elementBoundaryOffsets_new[sdN]+nElementBoundaries_subdomain_new[sdN];
898 valarray<int> elementBoundaryNumbering_new2old(elementBoundaries_subdomain_owned.size());
899 set<int>::iterator ebN_ownedp=elementBoundaries_subdomain_owned.begin();
900 for (
int ebN=0;ebN<int(elementBoundaries_subdomain_owned.size());ebN++)
902 elementBoundaryNumbering_new2old[ebN] = *ebN_ownedp++;
904 IS elementBoundaryNumberingIS_subdomain_new2old;
905 ISCreateGeneral(PROTEUS_COMM_WORLD,elementBoundaries_subdomain_owned.size(),&elementBoundaryNumbering_new2old[0],PETSC_COPY_VALUES,&elementBoundaryNumberingIS_subdomain_new2old);
906 IS elementBoundaryNumberingIS_global_new2old;
907 ISAllGather(elementBoundaryNumberingIS_subdomain_new2old,&elementBoundaryNumberingIS_global_new2old);
908 const PetscInt *elementBoundaryNumbering_global_new2old;
910 ISGetIndices(elementBoundaryNumberingIS_global_new2old,&elementBoundaryNumbering_global_new2old);
913 elementBoundaryNumbering_old2new_global[elementBoundaryNumbering_global_new2old[ebN]] = ebN;
927 map<NodeTuple<2>,
int> nodesEdgeMap_global;
928 set<int> edges_subdomain_owned;
933 const int nN0_global = nodeNumbering_global_old2new[nN0_global_old];
934 const int nN1_global = nodeNumbering_global_old2new[nN1_global_old];
936 nodes[0] = nN0_global;
937 nodes[1] = nN1_global;
939 nodesEdgeMap_global[et] = ig;
940 if (nodeOffsets_new[rank] <= et.
nodes[0] && et.
nodes[0] < nodeOffsets_new[rank+1])
941 edges_subdomain_owned.insert(ig);
944 valarray<int> nEdges_subdomain_new(size),
945 edgeOffsets_new(size+1);
947 for (
int sdN=0; sdN < size; sdN++)
949 nEdges_subdomain_new[sdN] = edges_subdomain_owned.size();
951 nEdges_subdomain_new[sdN] = 0;
953 valarray<int> nEdges_subdomain_new_send=nEdges_subdomain_new;
954 MPI_Allreduce(&nEdges_subdomain_new_send[0],&nEdges_subdomain_new[0],size,MPI_INT,MPI_SUM,PROTEUS_COMM_WORLD);
955 edgeOffsets_new[0] = 0;
956 for (
int sdN=0;sdN<size;sdN++)
957 edgeOffsets_new[sdN+1] = edgeOffsets_new[sdN]+nEdges_subdomain_new[sdN];
960 valarray<int> edgeNumbering_new2old(edges_subdomain_owned.size());
961 set<int>::iterator edges_ownedp = edges_subdomain_owned.begin();
962 for (
int i=0; i < int(edges_subdomain_owned.size());i++)
963 edgeNumbering_new2old[i] = *edges_ownedp++;
965 IS edgeNumberingIS_subdomain_new2old;
966 ISCreateGeneral(PROTEUS_COMM_WORLD,edges_subdomain_owned.size(),&edgeNumbering_new2old[0],PETSC_COPY_VALUES,&edgeNumberingIS_subdomain_new2old);
967 IS edgeNumberingIS_global_new2old;
968 ISAllGather(edgeNumberingIS_subdomain_new2old,&edgeNumberingIS_global_new2old);
969 const PetscInt *edgeNumbering_global_new2old;
970 valarray<int> edgeNumbering_old2new_global(mesh.
nEdges_global);
971 ISGetIndices(edgeNumberingIS_global_new2old,&edgeNumbering_global_new2old);
974 edgeNumbering_old2new_global[edgeNumbering_global_new2old[ig]] = ig;
980 valarray<int> edgeNodesArray_newNodesAndEdges(2*mesh.
nEdges_global);
981 map<NodeTuple<2>,
int > nodesEdgeMap_global_new;
986 const int nN0_global = nodeNumbering_global_old2new[nN0_global_old];
987 const int nN1_global = nodeNumbering_global_old2new[nN1_global_old];
989 const int edge_new = edgeNumbering_old2new_global[ig];
990 edgeNodesArray_newNodesAndEdges[edge_new*2+0] = nN0_global;
991 edgeNodesArray_newNodesAndEdges[edge_new*2+1] = nN1_global;
993 nodes[0] = nN0_global;
994 nodes[1] = nN1_global;
996 nodesEdgeMap_global_new[et] = edge_new;
1008 set<int> elements_overlap,nodes_overlap,elementBoundaries_overlap,edges_overlap;
1009 for (
int nN = nodeOffsets_new[rank]; nN < nodeOffsets_new[rank+1]; nN++)
1011 int nN_global_old = nodeNumbering_global_new2old[nN];
1017 int nN_neig_new = nodeNumbering_global_old2new[nN_neig_old];
1019 bool offproc = nN_neig_new >= nodeOffsets_new[rank+1] || nN_neig_new < nodeOffsets_new[rank];
1021 nodes_overlap.insert(nN_neig_new);
1028 int eN_star_new = elementNumbering_global_old2new[eN_star_old];
1029 bool offproc = eN_star_new >= elementOffsets_new[rank+1] || eN_star_new < elementOffsets_new[rank];
1031 elements_overlap.insert(eN_star_new);
1039 if (ebN_global < elementBoundaryOffsets_new[rank] || ebN_global >= elementBoundaryOffsets_new[rank+1])
1041 elementBoundaries_overlap.insert(ebN_global);
1049 const int nN0_global = nodeNumbering_global_old2new[nN0_global_old];
1050 const int nN1_global = nodeNumbering_global_old2new[nN1_global_old];
1051 bool foundEdge =
false;
1053 nodes[0] = nN0_global;
1054 nodes[1] = nN1_global;
1056 const int edge_global = nodesEdgeMap_global_new[et];
1057 if (edge_global < edgeOffsets_new[rank] || edge_global >= edgeOffsets_new[rank+1])
1058 edges_overlap.insert(edge_global);
1068 int overlap_remaining = nNodes_overlap -1;
1070 set<int> last_nodes_added2overlap = nodes_overlap;
1071 while (overlap_remaining > 0)
1073 set<int> new_nodes_overlap,new_elements_overlap,new_elementBoundaries_overlap,new_edges_overlap;
1074 set<int>::iterator nN_p = last_nodes_added2overlap.begin();
1075 while (nN_p != last_nodes_added2overlap.end())
1077 int nN_global_new = *nN_p;
1078 int nN_global_old = nodeNumbering_global_new2old[nN_global_new];
1083 int nN_neig_new = nodeNumbering_global_old2new[nN_neig_old];
1086 bool offproc = nN_neig_new >= nodeOffsets_new[rank+1] || nN_neig_new < nodeOffsets_new[rank];
1088 new_nodes_overlap.insert(nN_neig_new);
1094 set<int>::iterator nN_newp = last_nodes_added2overlap.begin();
1095 while (nN_newp != last_nodes_added2overlap.end())
1097 int nN_global_new = *nN_newp;
1098 int nN_global_old = nodeNumbering_global_new2old[nN_global_new];
1103 int eN_star_new = elementNumbering_global_old2new[eN_star_old];
1104 bool offproc = eN_star_new >= elementOffsets_new[rank+1] || eN_star_new < elementOffsets_new[rank];
1106 new_elements_overlap.insert(eN_star_new);
1115 if (ebN_global < elementBoundaryOffsets_new[rank] || ebN_global >= elementBoundaryOffsets_new[rank+1])
1116 new_elementBoundaries_overlap.insert(ebN_global);
1124 const int nN0_global = nodeNumbering_global_old2new[nN0_global_old];
1125 const int nN1_global = nodeNumbering_global_old2new[nN1_global_old];
1126 bool foundEdge =
false;
1128 nodes[0] = nN0_global;
1129 nodes[1] = nN1_global;
1131 const int edge_global = nodesEdgeMap_global_new[et];
1132 if (edge_global < edgeOffsets_new[rank] || edge_global >= edgeOffsets_new[rank+1])
1133 new_edges_overlap.insert(edge_global);
1138 last_nodes_added2overlap.clear();
1139 set_difference(new_nodes_overlap.begin(),new_nodes_overlap.end(),
1140 nodes_overlap.begin(),nodes_overlap.end(),
1141 insert_iterator<set<int> >(last_nodes_added2overlap,
1142 last_nodes_added2overlap.begin()));
1145 for (set<int>::iterator nN_addedp = new_nodes_overlap.begin();
1146 nN_addedp != new_nodes_overlap.end();
1149 nodes_overlap.insert(*nN_addedp);
1151 for (set<int>::iterator eN_addedp = new_elements_overlap.begin();
1152 eN_addedp != new_elements_overlap.end();
1155 elements_overlap.insert(*eN_addedp);
1157 for (set<int>::iterator ebN_addedp = new_elementBoundaries_overlap.begin();
1158 ebN_addedp != new_elementBoundaries_overlap.end();
1161 elementBoundaries_overlap.insert(*ebN_addedp);
1163 for (set<int>::iterator edge_addedp = new_edges_overlap.begin();
1164 edge_addedp != new_edges_overlap.end();
1167 edges_overlap.insert(*edge_addedp);
1171 overlap_remaining--;
1191 map<int,int> nodeNumbering_global2subdomain;
1192 map<int,int> elementBoundaryNumbering_global2subdomain;
1196 for (
int nN = 0; nN < nNodes_subdomain_new[rank]; nN++)
1198 int nN_global_new = nN + nodeOffsets_new[rank];
1199 int nN_global_old = nodeNumbering_global_new2old[nN_global_new];
1200 nodeNumbering_subdomain2global[nN] = nN_global_new;
1201 nodeNumbering_global2subdomain[nN_global_new] = nN;
1210 set<int>::iterator nN_p = nodes_overlap.begin();
1211 for (
int nN = nNodes_subdomain_new[rank]; nN < nNodes_subdomain_new[rank] + int(nodes_overlap.size()); nN++)
1213 int nN_global_new = *nN_p++;
1214 int nN_global_old = nodeNumbering_global_new2old[nN_global_new];
1215 nodeNumbering_subdomain2global[nN] = nN_global_new;
1216 nodeNumbering_global2subdomain[nN_global_new] = nN;
1225 for (
int eN = 0; eN < nElements_subdomain_new[rank]; eN++)
1227 int eN_global_new = elementOffsets_new[rank] + eN;
1228 int eN_global_old = elementNumbering_global_new2old[eN_global_new];
1229 elementNumbering_subdomain2global[eN] = eN_global_new;
1234 int nN_global_new = nodeNumbering_global_old2new[nN_global_old];
1235 int nN_subdomain = nodeNumbering_global2subdomain[nN_global_new];
1241 set<int>::iterator eN_p = elements_overlap.begin();
1242 for (
int eN = nElements_subdomain_new[rank]; eN < nElements_subdomain_new[rank] + int(elements_overlap.size()); eN++)
1244 int eN_global_new = *eN_p++;
1245 int eN_global_old = elementNumbering_global_new2old[eN_global_new];
1246 elementNumbering_subdomain2global[eN] = eN_global_new;
1251 int nN_global_new = nodeNumbering_global_old2new[nN_global_old];
1252 int nN_subdomain = nodeNumbering_global2subdomain[nN_global_new];
1258 for (
int ebN=0; ebN < nElementBoundaries_subdomain_new[rank]; ebN++)
1260 int ebN_global = ebN + elementBoundaryOffsets_new[rank];
1261 elementBoundaryNumbering_subdomain2global[ebN]=ebN_global;
1262 elementBoundaryNumbering_global2subdomain[ebN_global] = ebN;
1265 set<int>::iterator ebN_p = elementBoundaries_overlap.begin();
1266 for(
int ebN=nElementBoundaries_subdomain_new[rank];ebN < nElementBoundaries_subdomain_new[rank] + int(elementBoundaries_overlap.size()); ebN++)
1268 int ebN_global = *ebN_p++;
1269 elementBoundaryNumbering_subdomain2global[ebN] = ebN_global;
1270 elementBoundaryNumbering_global2subdomain[ebN_global] = ebN;
1275 for (
int eN=0;eN<nElements_subdomain_new[rank];eN++)
1277 int eN_global = eN+elementOffsets_new[rank];
1283 set<int>::iterator eN_p2 = elements_overlap.begin();
1284 for (
int eN = nElements_subdomain_new[rank]; eN < nElements_subdomain_new[rank] + int(elements_overlap.size()); eN++)
1286 int eN_global_new = *eN_p2++;
1296 for (
int i=0; i < nEdges_subdomain_new[rank]; i++)
1298 const int ig = i+edgeOffsets_new[rank];
1299 const int nN0_global = edgeNodesArray_newNodesAndEdges[ig*2+0];
1300 const int nN1_global = edgeNodesArray_newNodesAndEdges[ig*2+1];
1302 const int nN0_subdomain = nodeNumbering_global2subdomain[nN0_global];
1303 const int nN1_subdomain = nodeNumbering_global2subdomain[nN1_global];
1306 edgeNumbering_subdomain2global[i] = ig;
1309 set<int>::iterator edge_p = edges_overlap.begin();
1310 for (
int i=nEdges_subdomain_new[rank]; i < nEdges_subdomain_new[rank] + int(edges_overlap.size()); i++)
1312 const int ig =*edge_p++;
1313 const int nN0_global = edgeNodesArray_newNodesAndEdges[ig*2+0];
1314 const int nN1_global = edgeNodesArray_newNodesAndEdges[ig*2+1];
1316 const int nN0_subdomain = nodeNumbering_global2subdomain[nN0_global];
1317 const int nN1_subdomain = nodeNumbering_global2subdomain[nN1_global];
1320 edgeNumbering_subdomain2global[i] = ig;
1386 int eN_global_new = elementNumbering_subdomain2global[eN];
1387 int eN_global_old = elementNumbering_global_new2old[eN_global_new];
1409 for (
int sdN = 0; sdN < size+1; sdN++)
1440 ISRestoreIndices(nodeNumberingIS_global_old2new,&nodeNumbering_global_old2new);
1442 ISDestroy(&nodePartitioningIS_new);
1443 ISDestroy(&nodeNumberingIS_subdomain_old2new);
1444 ISDestroy(&nodeNumberingIS_global_old2new);
1446 ISRestoreIndices(elementNumberingIS_global_new2old,&elementNumbering_global_new2old);
1448 ISDestroy(&elementNumberingIS_subdomain_new2old);
1449 ISDestroy(&elementNumberingIS_global_new2old);
1451 ISRestoreIndices(elementBoundaryNumberingIS_global_new2old,&elementBoundaryNumbering_global_new2old);
1453 ISDestroy(&elementBoundaryNumberingIS_subdomain_new2old);
1454 ISDestroy(&elementBoundaryNumberingIS_global_new2old);
1456 ISRestoreIndices(edgeNumberingIS_global_new2old,&edgeNumbering_global_new2old);
1458 ISDestroy(&edgeNumberingIS_subdomain_new2old);
1459 ISDestroy(&edgeNumberingIS_global_new2old);
1466 using namespace std;
1467 PetscErrorCode ierr;
1468 PetscMPIInt size,rank;
1470 ierr = MPI_Comm_size(PROTEUS_COMM_WORLD,&size);CHKERRABORT(PROTEUS_COMM_WORLD, ierr);
1471 ierr = MPI_Comm_rank(PROTEUS_COMM_WORLD,&rank);CHKERRABORT(PROTEUS_COMM_WORLD, ierr);
1472 PetscLogStage partitioning_stage;
1473 PetscLogStageRegister(
"Mesh Partition",&partitioning_stage);
1474 PetscLogStagePush(partitioning_stage);
1501 bool failed =
false;
1502 const int simplexDim = 4;
1503 const int vertexDim = 3;
1505 std::string vertexFileName = std::string(filebase) +
".node" ;
1506 std::string elementFileName = std::string(filebase) +
".ele" ;
1507 std::string elementBoundaryFileName = std::string(filebase) +
".face" ;
1508 std::string edgeFileName = std::string(filebase) +
".edge" ;
1519 int read_elements_event;
1520 PetscLogEventRegister(
"Read eles",0,&read_elements_event);
1521 PetscLogEventBegin(read_elements_event,0,0,0,0);
1522 std::ifstream vertexFile(vertexFileName.c_str());
1523 if (!vertexFile.good())
1525 std::cerr<<
"cannot open Tetgen node file "
1526 <<vertexFileName<<std::endl;
1530 int hasVertexMarkers(0),hasVertexAttributes(0),nSpace(3),nNodes_global;
1532 vertexFile >>
eatcomments >> nNodes_global >> nSpace >> hasVertexAttributes >> hasVertexMarkers >>
eatline ;
1533 assert(nNodes_global > 0);
1534 assert(nSpace == 3);
1539 if (hasVertexAttributes > 0)
1541 std::cerr<<
"WARNING Tetgen nodes hasAttributes= "<<hasVertexAttributes
1542 <<
" > 0 will treat first value as integer id for boundary!!"<<std::endl;
1543 hasVertexMarkers = 1;
1549 valarray<int> nodeOffsets_old(size+1);
1550 nodeOffsets_old[0] = 0;
1551 for (
int sdN=0; sdN < size; sdN++)
1553 nodeOffsets_old[sdN+1] = nodeOffsets_old[sdN] +
1554 int(nNodes_global)/size + (int(nNodes_global)%size > sdN);
1556 int nNodes_subdomain_old = nodeOffsets_old[rank+1] - nodeOffsets_old[rank];
1563 std::ifstream elementFile(elementFileName.c_str());
1564 if (!elementFile.good())
1566 std::cerr<<
"cannot open Tetgen element file "
1567 <<elementFileName<<std::endl;
1572 int nNodesPerSimplex(simplexDim),hasElementMarkers = 0,nElements_global;
1573 elementFile >>
eatcomments >> nElements_global >> nNodesPerSimplex >> hasElementMarkers >>
eatline;
1574 assert(nElements_global > 0);
1575 assert(nNodesPerSimplex == simplexDim);
1577 vector<int> element_nodes_old(4);
1578 vector<set<int> > nodeStar(nNodes_subdomain_old);
1579 map<int,vector<int> > elements_old;
1580 for (
int ie = 0; ie < nElements_global; ie++)
1585 assert(0 <= ne && ne < nElements_global && elementFile.good());
1586 for (
int iv = 0; iv < simplexDim; iv++)
1590 assert(0 <= nv && nv < nNodes_global);
1591 element_nodes_old[iv] = nv;
1594 for (
int iv = 0; iv < simplexDim; iv++)
1597 int nN_star = element_nodes_old[iv];
1598 bool inSubdomain=
false;
1599 if (nN_star >= nodeOffsets_old[rank] && nN_star < nodeOffsets_old[rank+1])
1603 for (
int jv = 0; jv < simplexDim; jv++)
1607 int nN_star_subdomain = nN_star-nodeOffsets_old[rank];
1608 nodeStar[nN_star_subdomain].insert(element_nodes_old[jv]);
1613 elements_old[ie] = element_nodes_old;
1617 elementFile.close();
1618 PetscLogEventEnd(read_elements_event,0,0,0,0);
1619 int repartition_nodes_event;
1620 PetscLogEventRegister(
"Repart nodes",0,&repartition_nodes_event);
1621 PetscLogEventBegin(repartition_nodes_event,0,0,0,0);
1624 valarray<int> nodeStarOffsets(nNodes_subdomain_old+1);
1625 nodeStarOffsets[0] = 0;
1626 for (
int nN=1;nN<nNodes_subdomain_old+1;nN++)
1627 nodeStarOffsets[nN] = nodeStarOffsets[nN-1] + nodeStar[nN-1].size();
1628 valarray<int> nodeStarArray(nodeStarOffsets[nNodes_subdomain_old]);
1629 for (
int nN=0,offset=0;nN<nNodes_subdomain_old;nN++)
1630 for (set<int>::iterator nN_star=nodeStar[nN].begin();nN_star!=nodeStar[nN].end();nN_star++,offset++)
1631 nodeStarArray[offset] = *nN_star;
1633 int max_nNodeNeighbors_node=0;
1634 for (
int nN=0;nN<nNodes_subdomain_old;nN++)
1635 max_nNodeNeighbors_node=
max(max_nNodeNeighbors_node,nodeStarOffsets[nN+1]-nodeStarOffsets[nN]);
1637 PetscBool isInitialized;
1638 PetscInitialized(&isInitialized);
1639 PetscInt *nodeNeighborsOffsets_subdomain,*nodeNeighbors_subdomain,*weights_subdomain,*vertex_weights_subdomain;
1640 PetscReal *partition_weights;
1641 PetscMalloc(
sizeof(PetscInt)*(nNodes_subdomain_old+1),&nodeNeighborsOffsets_subdomain);
1642 PetscMalloc(
sizeof(PetscInt)*(nNodes_subdomain_old*max_nNodeNeighbors_node),&nodeNeighbors_subdomain);
1643 PetscMalloc(
sizeof(PetscInt)*(nNodes_subdomain_old*max_nNodeNeighbors_node),&weights_subdomain);
1644 PetscMalloc(
sizeof(PetscInt)*(nNodes_subdomain_old),&vertex_weights_subdomain);
1645 PetscMalloc(
sizeof(PetscReal)*(size),&partition_weights);
1646 for (
int sd=0;sd<size;sd++)
1647 partition_weights[sd] = 1.0/
double(size);
1648 nodeNeighborsOffsets_subdomain[0] = 0;
1650 for (
int nN = 0,offset=0; nN < nNodes_subdomain_old; nN++)
1652 for (
int offset_subdomain = nodeStarOffsets[nN];
1653 offset_subdomain < nodeStarOffsets[nN+1];
1656 nodeNeighbors_subdomain[offset++] = nodeStarArray[offset_subdomain];
1658 nodeNeighborsOffsets_subdomain[nN+1]=offset;
1659 sort(&nodeNeighbors_subdomain[nodeNeighborsOffsets_subdomain[nN]],&nodeNeighbors_subdomain[nodeNeighborsOffsets_subdomain[nN+1]]);
1661 int weight= (nodeNeighborsOffsets_subdomain[nN+1] - nodeNeighborsOffsets_subdomain[nN]);
1662 vertex_weights_subdomain[nN] = weight;
1663 for (
int k=nodeNeighborsOffsets_subdomain[nN];k<nodeNeighborsOffsets_subdomain[nN+1];k++)
1664 weights_subdomain[k] = weight;
1670 int nNodes_subdomain_max=0;
1671 MPI_Allreduce(&nNodes_subdomain_old,
1672 &nNodes_subdomain_max,
1676 PROTEUS_COMM_WORLD);
1678 std::cout<<
"Max nNodes_subdomain "<<nNodes_subdomain_max<<
" nNodes_global "<<nNodes_global<<std::endl;
1679 ierr = MatCreateMPIAdj(PROTEUS_COMM_WORLD,
1680 nNodes_subdomain_old,
1682 nodeNeighborsOffsets_subdomain,
1683 nodeNeighbors_subdomain,
1685 &petscAdjacency);CHKERRABORT(PROTEUS_COMM_WORLD, ierr);
1687 const double max_rss_gb(3.75);
1688 ierr =
enforceMemoryLimit(PROTEUS_COMM_WORLD, rank, max_rss_gb,
"Done allocating MPIAdj");CHKERRABORT(PROTEUS_COMM_WORLD, ierr);
1689 MatPartitioning petscPartition;
1690 ierr = MatPartitioningCreate(PROTEUS_COMM_WORLD,&petscPartition);CHKERRABORT(PROTEUS_COMM_WORLD, ierr);
1691 ierr = MatPartitioningSetAdjacency(petscPartition,petscAdjacency);CHKERRABORT(PROTEUS_COMM_WORLD, ierr);
1692 ierr = MatPartitioningSetFromOptions(petscPartition);CHKERRABORT(PROTEUS_COMM_WORLD, ierr);
1693 ierr = MatPartitioningSetVertexWeights(petscPartition,vertex_weights_subdomain);CHKERRABORT(PROTEUS_COMM_WORLD, ierr);
1694 ierr = MatPartitioningSetPartitionWeights(petscPartition,partition_weights);CHKERRABORT(PROTEUS_COMM_WORLD, ierr);
1696 IS nodePartitioningIS_new;
1697 ierr = MatPartitioningApply(petscPartition,&nodePartitioningIS_new);CHKERRABORT(PROTEUS_COMM_WORLD, ierr);
1698 ierr = MatPartitioningDestroy(&petscPartition);CHKERRABORT(PROTEUS_COMM_WORLD, ierr);
1699 ierr =
enforceMemoryLimit(PROTEUS_COMM_WORLD, rank, max_rss_gb,
"Done applying partition");CHKERRABORT(PROTEUS_COMM_WORLD, ierr);
1702 valarray<int> nNodes_subdomain_new(size);
1703 ISPartitioningCount(nodePartitioningIS_new,size,&nNodes_subdomain_new[0]);
1706 valarray<int> nodeOffsets_new(size+1);
1707 nodeOffsets_new[0] = 0;
1708 for (
int sdN = 0; sdN < size; sdN++)
1710 nodeOffsets_new[sdN+1] = nodeOffsets_new[sdN] + nNodes_subdomain_new[sdN];
1714 IS nodeNumberingIS_subdomain_old2new;
1715 ISPartitioningToNumbering(nodePartitioningIS_new,&nodeNumberingIS_subdomain_old2new);
1722 MPI_Info info = MPI_INFO_NULL;
1723 hid_t plist_id = H5Pcreate(H5P_FILE_ACCESS);
1724 H5Pset_fapl_mpio(plist_id, PROTEUS_COMM_WORLD, info);
1729 const char* H5FILE_NAME(
"mappings.h5");
1730 hid_t file_id = H5Fcreate(H5FILE_NAME, H5F_ACC_TRUNC, H5P_DEFAULT, plist_id);
1738 dimsf[0] = nNodes_global;
1740 hid_t filespace = H5Screate_simple(
RANK, dimsf, NULL);
1745 hid_t dset_id = H5Dcreate(file_id,
"nodeNumbering_old2new", H5T_NATIVE_INT, filespace,
1746 H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
1747 H5Sclose(filespace);
1755 count[0] = nNodes_subdomain_old;
1756 offset[0] = nodeOffsets_old[rank];
1757 hid_t memspace = H5Screate_simple(
RANK, count, NULL);
1762 filespace = H5Dget_space(dset_id);
1763 H5Sselect_hyperslab(filespace, H5S_SELECT_SET, offset, NULL, count, NULL);
1772 const PetscInt* data;
1773 ISGetIndices(nodeNumberingIS_subdomain_old2new, &data);
1778 plist_id = H5Pcreate(H5P_DATASET_XFER);
1779 H5Pset_dxpl_mpio(plist_id, H5FD_MPIO_COLLECTIVE);
1781 herr_t status = H5Dwrite(dset_id, H5T_NATIVE_INT, memspace, filespace,
1784 ISRestoreIndices(nodeNumberingIS_subdomain_old2new, &data);
1795 IS nodeNumberingIS_global_old2new;
1796 ISAllGather(nodeNumberingIS_subdomain_old2new,&nodeNumberingIS_global_old2new);
1797 const PetscInt * nodeNumbering_global_old2new;
1798 ISGetIndices(nodeNumberingIS_global_old2new,&nodeNumbering_global_old2new);
1806 int dset_data[nNodes_global];
1812 dataset_id = H5Dopen2(file_id,
"/nodeNumbering_old2new", H5P_DEFAULT);
1814 status = H5Dread(dataset_id, H5T_NATIVE_INT, H5S_ALL, H5S_ALL, H5P_DEFAULT,
1818 status = H5Dclose(dataset_id);
1820 for (
int i=0;i<nNodes_global;i++)
1821 assert(nodeNumbering_global_old2new[i] == dset_data[i]);
1822 std::cout<<
"==================out of core old2new is correct!===================="<<std::endl;
1834 PetscLogEventEnd(repartition_nodes_event,0,0,0,0);
1835 int receive_element_mask_event;
1836 PetscLogEventRegister(
"Recv. ele mask",0,&receive_element_mask_event);
1837 PetscLogEventBegin(receive_element_mask_event,0,0,0,0);
1842 PetscLogEventEnd(receive_element_mask_event,0,0,0,0);
1843 ierr =
enforceMemoryLimit(PROTEUS_COMM_WORLD, rank, max_rss_gb,
"Done with masks");CHKERRABORT(PROTEUS_COMM_WORLD, ierr);
1844 int build_subdomains_reread_elements_event;
1845 PetscLogEventRegister(
"Reread eles",0,&build_subdomains_reread_elements_event);
1846 PetscLogEventBegin(build_subdomains_reread_elements_event,0,0,0,0);
1853 std::ifstream elementFile2(elementFileName.c_str());
1855 if (!elementFile2.good())
1857 std::cerr<<
"cannot open Tetgen elements file"
1858 <<elementFileName<<std::endl;
1862 elementFile2 >>
eatcomments >> nElements_global >> nNodesPerSimplex >> hasElementMarkers >>
eatline;
1863 assert(nElements_global > 0);
1864 assert(nNodesPerSimplex == simplexDim);
1865 set<int> elements_subdomain_owned;
1866 vector<int> element_nodes_new(4);
1867 int element_nodes_new_array[4];
1868 vector<set<int> > nodeElementsStar(nNodes_subdomain_new[rank]);
1869 vector<set<int> > nodeStarNew(nNodes_subdomain_new[rank]);
1870 map<int,vector<int> > elementNodesArrayMap;
1871 map<int,long int> elementMaterialTypesMap;
1873 map<NodeTuple<2>,set<pair<int,int> > > edgeElementsMap;
1875 for (
int ie = 0; ie < nElements_global; ie++)
1877 int ne, nv, elementId(0);
1878 long double elementId_double;
1881 assert(0 <= ne && ne < nElements_global && elementFile.good());
1882 for (
int iv = 0; iv < simplexDim; iv++)
1884 elementFile2 >> nv ;
1886 assert(0 <= nv && nv < nNodes_global);
1887 element_nodes_old[iv] = nv;
1888 element_nodes_new[iv] = nodeNumbering_global_old2new[nv];
1889 element_nodes_new_array[iv] = element_nodes_new[iv];
1892 for (
int iv = 0; iv < simplexDim; iv++)
1894 int nN_star_new = element_nodes_new[iv];
1895 bool inSubdomain=
false;
1896 if (nN_star_new >= nodeOffsets_new[rank] && nN_star_new < nodeOffsets_new[rank+1])
1900 for (
int ebN=0;ebN < 4 ; ebN++)
1902 int nodes[3] = { element_nodes_new[(ebN+1) % 4],
1903 element_nodes_new[(ebN+2) % 4],
1904 element_nodes_new[(ebN+3) % 4]};
1906 if(elementBoundaryElementsMap.find(nodeTuple) != elementBoundaryElementsMap.end())
1908 if (elementBoundaryElementsMap[nodeTuple].right == -1 && ne != elementBoundaryElementsMap[nodeTuple].left)
1910 elementBoundaryElementsMap[nodeTuple].right=ne;
1911 elementBoundaryElementsMap[nodeTuple].right_ebN_element=ebN;
1920 for (
int nNL=0,edN=0;nNL < 4 ; nNL++)
1921 for(
int nNR=nNL+1;nNR < 4;nNR++,edN++)
1923 int nodes[2] = { element_nodes_new[nNL],
1924 element_nodes_new[nNR]};
1926 edgeElementsMap[nodeTuple].insert(pair<int,int>(ne,edN));
1929 int nN_star_new_subdomain = nN_star_new - nodeOffsets_new[rank];
1930 nodeElementsStar[nN_star_new_subdomain].insert(ne);
1931 for (
int jv = 0; jv < simplexDim; jv++)
1935 int nN_point_new = element_nodes_new[jv];
1936 nodeStarNew[nN_star_new_subdomain].insert(nN_point_new);
1942 elementNodesArrayMap[ne] = element_nodes_new;
1945 if (elementNodesArrayMap.find(ne) != elementNodesArrayMap.end())
1947 if (nodeTuple.
nodes[1] >= nodeOffsets_new[rank] && nodeTuple.
nodes[1] < nodeOffsets_new[rank+1])
1948 elements_subdomain_owned.insert(ne);
1949 if (hasElementMarkers > 0)
1951 elementFile2 >> elementId_double;
1952 elementId =
static_cast<long int>(elementId_double);
1953 elementMaterialTypesMap[ne] = elementId;
1958 elementFile2.close();
1959 int nElements_owned_subdomain(elements_subdomain_owned.size()),
1960 nElements_owned_new=0;
1961 MPI_Allreduce(&nElements_owned_subdomain,&nElements_owned_new,1,MPI_INT,MPI_SUM,PROTEUS_COMM_WORLD);
1962 assert(nElements_owned_new == nElements_global);
1963 PetscLogEventEnd(build_subdomains_reread_elements_event,0,0,0,0);
1964 int build_subdomains_send_marked_elements_event;
1965 PetscLogEventRegister(
"Mark/send eles",0,&build_subdomains_send_marked_elements_event);
1966 PetscLogEventBegin(build_subdomains_send_marked_elements_event,0,0,0,0);
1967 ierr =
enforceMemoryLimit(PROTEUS_COMM_WORLD, rank, max_rss_gb,
"Done marking elements");CHKERRABORT(PROTEUS_COMM_WORLD, ierr);
1972 valarray<int> nodeElementOffsets(nNodes_subdomain_new[rank]+1);
1973 nodeElementOffsets[0] = 0;
1974 for (
int nN = 0; nN < nNodes_subdomain_new[rank]; nN++)
1975 nodeElementOffsets[nN+1] = nodeElementOffsets[nN]+nodeElementsStar[nN].size();
1976 valarray<int> nodeElementsArray(nodeElementOffsets[nNodes_subdomain_new[rank]]);
1977 for (
int nN=0,offset=0; nN < nNodes_subdomain_new[rank]; nN++)
1979 for (set<int>::iterator eN_star = nodeElementsStar[nN].begin(); eN_star != nodeElementsStar[nN].end();
1982 nodeElementsArray[offset] = *eN_star;
1986 valarray<int> nodeStarOffsetsNew(nNodes_subdomain_new[rank]+1);
1987 nodeStarOffsetsNew[0] = 0;
1988 for (
int nN=1;nN<nNodes_subdomain_new[rank]+1;nN++)
1989 nodeStarOffsetsNew[nN] = nodeStarOffsetsNew[nN-1] + nodeStarNew[nN-1].size();
1990 valarray<int> nodeStarArrayNew(nodeStarOffsetsNew[nNodes_subdomain_new[rank]]);
1991 for (
int nN=0,offset=0;nN<nNodes_subdomain_new[rank];nN++)
1993 for (set<int>::iterator nN_star=nodeStarNew[nN].begin();nN_star!=nodeStarNew[nN].end();nN_star++,offset++)
1995 nodeStarArrayNew[offset] = *nN_star;
1998 PetscLogEventEnd(build_subdomains_send_marked_elements_event,0,0,0,0);
1999 int build_subdomains_global_numbering_elements_event;
2000 PetscLogEventRegister(
"Global ele nmbr",0,&build_subdomains_global_numbering_elements_event);
2001 PetscLogEventBegin(build_subdomains_global_numbering_elements_event,0,0,0,0);
2005 valarray<int> nElements_subdomain_new(size),
2006 elementOffsets_new(size+1);
2007 for (
int sdN = 0; sdN < size; sdN++)
2010 nElements_subdomain_new[sdN] = int(elements_subdomain_owned.size());
2012 nElements_subdomain_new[sdN] = 0;
2014 valarray<int> nElements_subdomain_new_send = nElements_subdomain_new;
2015 MPI_Allreduce(&nElements_subdomain_new_send[0],&nElements_subdomain_new[0],size,MPI_INT,MPI_SUM,PROTEUS_COMM_WORLD);
2017 elementOffsets_new[0] = 0;
2018 for (
int sdN = 0; sdN < size; sdN++)
2019 elementOffsets_new[sdN+1] = elementOffsets_new[sdN] + nElements_subdomain_new[sdN];
2021 valarray<int> elementNumbering_subdomain_new2old(elements_subdomain_owned.size());
2022 set<int>::iterator eN_ownedp = elements_subdomain_owned.begin();
2023 for (
int eN = 0; eN < int(elements_subdomain_owned.size()); eN++,eN_ownedp++)
2025 elementNumbering_subdomain_new2old[eN] = *eN_ownedp;
2028 IS elementNumberingIS_subdomain_new2old;
2029 ISCreateGeneral(PROTEUS_COMM_WORLD,elements_subdomain_owned.size(),&elementNumbering_subdomain_new2old[0],PETSC_COPY_VALUES,
2030 &elementNumberingIS_subdomain_new2old);
2031 IS elementNumberingIS_global_new2old;
2032 ISAllGather(elementNumberingIS_subdomain_new2old,&elementNumberingIS_global_new2old);
2034 const PetscInt *elementNumbering_global_new2old;
2035 ISGetIndices(elementNumberingIS_global_new2old,&elementNumbering_global_new2old);
2037 valarray<int> elementNumbering_global_old2new(nElements_global);
2038 for (
int eN = 0; eN < nElements_global; eN++)
2040 elementNumbering_global_old2new[elementNumbering_global_new2old[eN]] = eN;
2042 PetscLogEventEnd(build_subdomains_global_numbering_elements_event,0,0,0,0);
2043 ierr =
enforceMemoryLimit(PROTEUS_COMM_WORLD, rank, max_rss_gb,
"Done allocating element numbering new2old/old2new");CHKERRABORT(PROTEUS_COMM_WORLD, ierr);
2044 int build_subdomains_faces_event;
2045 PetscLogEventRegister(
"Subd faces",0,&build_subdomains_faces_event);
2046 PetscLogEventBegin(build_subdomains_faces_event,0,0,0,0);
2054 std::ifstream elementBoundaryFile(elementBoundaryFileName.c_str());
2056 if (!elementBoundaryFile.good())
2058 std::cerr<<
"cannot open Tetgen face file "
2059 <<elementBoundaryFileName<<std::endl;
2064 bool hasElementBoundaryMarkers =
false;
2065 int nElementBoundaries_global;
2066 int ihasElementBoundaryMarkers(0);
2067 elementBoundaryFile >>
eatcomments >> nElementBoundaries_global >> ihasElementBoundaryMarkers >>
eatline ;
2068 assert(nElementBoundaries_global > 0);
2069 if (ihasElementBoundaryMarkers > 0)
2071 hasElementBoundaryMarkers =
true;
2075 set<int> elementBoundaries_subdomain_owned;
2076 vector<set<int> > nodeElementBoundariesStar(nNodes_subdomain_new[rank]);
2077 map<int,int> elementBoundaryMaterialTypesMap;
2078 map<int,vector<int> > elementBoundariesMap;
2079 set<int> supportedElementBoundaries;
2080 for (
int ieb = 0; ieb < nElementBoundaries_global; ieb++)
2082 int neb,nn0,nn1,nn2;
int ebId(0);
2083 elementBoundaryFile >>
eatcomments >> neb >> nn0 >> nn1 >> nn2;
2084 if (ihasElementBoundaryMarkers > 0)
2085 elementBoundaryFile >> ebId;
2090 assert(0 <= neb && neb < nElementBoundaries_global && elementBoundaryFile.good());
2093 int nn0_new = nodeNumbering_global_old2new[nn0];
2094 if (nn0_new >= nodeOffsets_new[rank] && nn0_new < nodeOffsets_new[rank+1])
2096 nodeElementBoundariesStar[nn0_new-nodeOffsets_new[rank]].insert(neb);
2097 supportedElementBoundaries.insert(neb);
2099 int nn1_new = nodeNumbering_global_old2new[nn1];
2100 if (nn1_new >= nodeOffsets_new[rank] && nn1_new < nodeOffsets_new[rank+1])
2102 nodeElementBoundariesStar[nn1_new-nodeOffsets_new[rank]].insert(neb);
2103 supportedElementBoundaries.insert(neb);
2105 int nn2_new = nodeNumbering_global_old2new[nn2];
2106 if (nn2_new >= nodeOffsets_new[rank] && nn2_new < nodeOffsets_new[rank+1])
2108 nodeElementBoundariesStar[nn2_new-nodeOffsets_new[rank]].insert(neb);
2109 supportedElementBoundaries.insert(neb);
2111 int nodes[3] = {nn0_new,nn1_new,nn2_new};
2113 elementBoundaryFile >>
eatline;
2114 if (elementBoundaryElementsMap.find(nodeTuple) != elementBoundaryElementsMap.end())
2116 if (nodeTuple.
nodes[1] >= nodeOffsets_new[rank] && nodeTuple.
nodes[1] < nodeOffsets_new[rank+1])
2117 elementBoundaries_subdomain_owned.insert(neb);
2118 if (ihasElementBoundaryMarkers > 0)
2119 elementBoundaryMaterialTypesMap[neb]=ebId;
2120 int eN_left = elementNumbering_global_old2new[elementBoundaryElementsMap[nodeTuple].left];
2121 if (elementBoundariesMap.find(eN_left) != elementBoundariesMap.end())
2123 elementBoundariesMap[eN_left][elementBoundaryElementsMap[nodeTuple].left_ebN_element] = neb;
2128 vector<int> elementBoundaries_element(4,-1);
2129 elementBoundariesMap[eN_left] = elementBoundaries_element;
2131 elementBoundariesMap[eN_left][elementBoundaryElementsMap[nodeTuple].left_ebN_element] = neb;
2133 if (elementBoundaryElementsMap[nodeTuple].right >= 0)
2135 int eN_right = elementNumbering_global_old2new[elementBoundaryElementsMap[nodeTuple].right];
2136 if (elementBoundariesMap.find(eN_right) != elementBoundariesMap.end())
2138 elementBoundariesMap[eN_right][elementBoundaryElementsMap[nodeTuple].right_ebN_element] = neb;
2143 vector<int> elementBoundaries_element(4,-1);
2144 elementBoundariesMap[eN_right] = elementBoundaries_element;
2146 elementBoundariesMap[eN_right][elementBoundaryElementsMap[nodeTuple].right_ebN_element] = neb;
2152 elementBoundaryFile.close();
2153 ierr =
enforceMemoryLimit(PROTEUS_COMM_WORLD, rank, max_rss_gb,
"Done reading element boundaries");CHKERRABORT(PROTEUS_COMM_WORLD, ierr);
2154 int nElementBoundaries_owned_subdomain=elementBoundaries_subdomain_owned.size(),
2155 nElementBoundaries_owned_new=0;
2156 MPI_Allreduce(&nElementBoundaries_owned_subdomain,&nElementBoundaries_owned_new,1,MPI_INT,MPI_SUM,PROTEUS_COMM_WORLD);
2157 assert(nElementBoundaries_owned_new == nElementBoundaries_global);
2160 for (map<
int,vector<int> >::iterator elementBoundariesp=elementBoundariesMap.begin();
2161 elementBoundariesp!=elementBoundariesMap.end();
2162 elementBoundariesp++)
2165 for (
int iv=0;iv<4;iv++)
2168 int nN_global = elementNodesArrayMap[elementNumbering_global_new2old[elementBoundariesp->first]][iv];
2169 if (nN_global >= nodeOffsets_new[rank] && nN_global < nodeOffsets_new[rank+1])
2172 for(
int eb=0;eb<4;eb++)
2174 nodeElementBoundariesStar[nN_global-nodeOffsets_new[rank]].insert(elementBoundariesp->second[eb]);
2180 valarray<int> nodeElementBoundaryOffsets(nNodes_subdomain_new[rank]+1);
2181 nodeElementBoundaryOffsets[0] = 0;
2182 for (
int nN = 0; nN < nNodes_subdomain_new[rank]; nN++)
2183 nodeElementBoundaryOffsets[nN+1] = nodeElementBoundaryOffsets[nN]+nodeElementBoundariesStar[nN].size();
2184 valarray<int> nodeElementBoundariesArray(nodeElementBoundaryOffsets[nNodes_subdomain_new[rank]]);
2185 for (
int nN=0,offset=0; nN < nNodes_subdomain_new[rank]; nN++)
2187 for (set<int>::iterator ebN_star = nodeElementBoundariesStar[nN].begin(); ebN_star != nodeElementBoundariesStar[nN].end();
2188 ebN_star++,offset++)
2190 nodeElementBoundariesArray[offset] = *ebN_star;
2195 valarray<int> nElementBoundaries_subdomain_new(size),
2196 elementBoundaryOffsets_new(size+1);
2197 for (
int sdN=0;sdN<size;sdN++)
2199 nElementBoundaries_subdomain_new[sdN] = elementBoundaries_subdomain_owned.size();
2201 nElementBoundaries_subdomain_new[sdN] = 0;
2202 valarray<int> nElementBoundaries_subdomain_new_send=nElementBoundaries_subdomain_new;
2203 MPI_Allreduce(&nElementBoundaries_subdomain_new_send[0],&nElementBoundaries_subdomain_new[0],size,MPI_INT,MPI_SUM,PROTEUS_COMM_WORLD);
2204 elementBoundaryOffsets_new[0] = 0;
2205 for (
int sdN=0;sdN<size;sdN++)
2206 elementBoundaryOffsets_new[sdN+1] = elementBoundaryOffsets_new[sdN]+nElementBoundaries_subdomain_new[sdN];
2212 valarray<int> elementBoundaryNumbering_new2old(elementBoundaries_subdomain_owned.size());
2213 set<int>::iterator ebN_ownedp=elementBoundaries_subdomain_owned.begin();
2214 for (
int ebN=0;ebN<int(elementBoundaries_subdomain_owned.size());ebN++)
2216 elementBoundaryNumbering_new2old[ebN] = *ebN_ownedp++;
2218 IS elementBoundaryNumberingIS_subdomain_new2old;
2219 ISCreateGeneral(PROTEUS_COMM_WORLD,elementBoundaries_subdomain_owned.size(),&elementBoundaryNumbering_new2old[0],PETSC_COPY_VALUES,&elementBoundaryNumberingIS_subdomain_new2old);
2220 IS elementBoundaryNumberingIS_global_new2old;
2221 ISAllGather(elementBoundaryNumberingIS_subdomain_new2old,&elementBoundaryNumberingIS_global_new2old);
2222 const PetscInt *elementBoundaryNumbering_global_new2old;
2224 ISGetIndices(elementBoundaryNumberingIS_global_new2old,&elementBoundaryNumbering_global_new2old);
2225 ierr =
enforceMemoryLimit(PROTEUS_COMM_WORLD, rank, max_rss_gb,
"Allocating elementBoudnary old2new/new2old");CHKERRABORT(PROTEUS_COMM_WORLD, ierr);
2228 elementBoundaryNumbering_global_old2new[elementBoundaryNumbering_global_new2old[ebN]] = ebN;
2230 ISRestoreIndices(elementBoundaryNumberingIS_global_new2old,&elementBoundaryNumbering_global_new2old);
2231 ISDestroy(&elementBoundaryNumberingIS_subdomain_new2old);
2232 ISDestroy(&elementBoundaryNumberingIS_global_new2old);
2233 PetscLogEventEnd(build_subdomains_faces_event,0,0,0,0);
2234 ierr =
enforceMemoryLimit(PROTEUS_COMM_WORLD, rank, max_rss_gb,
"Done allocating elementBoudnary old2new/new2old");CHKERRABORT(PROTEUS_COMM_WORLD, ierr);
2235 int build_subdomains_edges_event;
2236 PetscLogEventRegister(
"Subd edges",0,&build_subdomains_edges_event);
2237 PetscLogEventBegin(build_subdomains_edges_event,0,0,0,0);
2242 std::ifstream edgeFile(edgeFileName.c_str());
2244 if (!edgeFile.good())
2246 std::cerr<<
"cannot open Tetgen edge file"
2247 <<edgeFileName<<std::endl;
2252 bool hasEdgeMarkers =
false;
2254 int ihasEdgeMarkers(0);
2256 assert(nEdges_global > 0);
2257 if (ihasEdgeMarkers > 0)
2259 hasEdgeMarkers =
true;
2262 set<int> edges_subdomain_owned;
2263 vector<set<int> > nodeEdgesStar(nNodes_subdomain_new[rank]);
2264 map<int,int> edgeMaterialTypesMap;
2265 map<int,vector<int> > elementEdgesMap;
2266 map<int,pair<int,int> > edgeNodesMap;
2267 set<int> supportedEdges;
2268 for (
int ied = 0; ied < nEdges_global; ied++)
2270 int ned,nn0,nn1;
int edId(0);
2272 if (ihasEdgeMarkers > 0)
2277 assert(0 <= ned && ned < nEdges_global && edgeFile.good());
2278 int nn0_new = nodeNumbering_global_old2new[nn0];
2279 if (nn0_new >= nodeOffsets_new[rank] && nn0_new < nodeOffsets_new[rank+1])
2281 nodeEdgesStar.at(nn0_new-nodeOffsets_new[rank]).insert(ned);
2282 supportedEdges.insert(ned);
2284 int nn1_new = nodeNumbering_global_old2new[nn1];
2285 if (nn1_new >= nodeOffsets_new[rank] && nn1_new < nodeOffsets_new[rank+1])
2287 nodeEdgesStar.at(nn1_new-nodeOffsets_new[rank]).insert(ned);
2288 supportedEdges.insert(ned);
2290 int nodes[2] = {nn0_new,nn1_new};
2293 if (edgeElementsMap.find(nodeTuple) != edgeElementsMap.end())
2295 if (nodeTuple.
nodes[0] >= nodeOffsets_new[rank] && nodeTuple.
nodes[0] < nodeOffsets_new[rank+1])
2296 edges_subdomain_owned.insert(ned);
2298 edgeNodesMap[ned].first = nodeTuple.
nodes[0];
2299 edgeNodesMap[ned].second = nodeTuple.
nodes[1];
2300 if (ihasEdgeMarkers > 0)
2301 edgeMaterialTypesMap[ned]=edId;
2302 for (set<pair<int,int> >::iterator elementp=edgeElementsMap[nodeTuple].begin();
2303 elementp != edgeElementsMap[nodeTuple].end();
2306 int eN = elementNumbering_global_old2new[elementp->first];
2307 if (elementEdgesMap.find(eN) != elementEdgesMap.end())
2309 elementEdgesMap[eN][elementp->second] = ned;
2313 std::vector<int>
init(6,-1);
2314 elementEdgesMap[eN] =
init;
2315 elementEdgesMap[eN][elementp->second] = ned;
2321 int nEdges_owned_subdomain=edges_subdomain_owned.size(),
2323 MPI_Allreduce(&nEdges_owned_subdomain,&nEdges_owned_new,1,MPI_INT,MPI_SUM,PROTEUS_COMM_WORLD);
2324 assert(nEdges_owned_new == nEdges_global);
2325 ierr =
enforceMemoryLimit(PROTEUS_COMM_WORLD, rank, max_rss_gb,
"Done reading edges");CHKERRABORT(PROTEUS_COMM_WORLD, ierr);
2331 for (map<
int,vector<int> >::iterator edgesp=elementEdgesMap.begin();
2332 edgesp!=elementEdgesMap.end();
2336 for (
int iv=0;iv<4;iv++)
2339 int nN_global = elementNodesArrayMap[elementNumbering_global_new2old[edgesp->first]][iv];
2340 if (nN_global >= nodeOffsets_new[rank] && nN_global < nodeOffsets_new[rank+1])
2343 for(
int ed=0;ed<6;ed++)
2345 nodeEdgesStar.at(nN_global-nodeOffsets_new[rank]).insert(edgesp->second[ed]);
2351 valarray<int> nodeEdgeOffsets(nNodes_subdomain_new[rank]+1);
2352 nodeEdgeOffsets[0] = 0;
2353 for (
int nN = 0; nN < nNodes_subdomain_new[rank]; nN++)
2354 nodeEdgeOffsets[nN+1] = nodeEdgeOffsets[nN]+nodeEdgesStar.at(nN).size();
2355 valarray<int> nodeEdgesArray(nodeEdgeOffsets[nNodes_subdomain_new[rank]]);
2356 for (
int nN=0,offset=0; nN < nNodes_subdomain_new[rank]; nN++)
2358 for (set<int>::iterator edN_star = nodeEdgesStar.at(nN).begin();
2359 edN_star != nodeEdgesStar.at(nN).end();
2360 edN_star++,offset++)
2362 nodeEdgesArray[offset] = *edN_star;
2367 valarray<int> nEdges_subdomain_new(size),
2368 edgeOffsets_new(size+1);
2369 for (
int sdN=0;sdN<size;sdN++)
2371 nEdges_subdomain_new[sdN] = edges_subdomain_owned.size();
2373 nEdges_subdomain_new[sdN] = 0;
2374 valarray<int> nEdges_subdomain_new_send=nEdges_subdomain_new;
2375 MPI_Allreduce(&nEdges_subdomain_new_send[0],&nEdges_subdomain_new[0],size,MPI_INT,MPI_SUM,PROTEUS_COMM_WORLD);
2379 edgeOffsets_new[0] = 0;
2380 for (
int sdN=0;sdN<size;sdN++)
2381 edgeOffsets_new[sdN+1] = edgeOffsets_new[sdN]+nEdges_subdomain_new[sdN];
2387 valarray<int> edgeNumbering_new2old(edges_subdomain_owned.size());
2388 set<int>::iterator edN_ownedp=edges_subdomain_owned.begin();
2389 for (
int edN=0;edN<int(edges_subdomain_owned.size());edN++,edN_ownedp++)
2391 edgeNumbering_new2old[edN] = *edN_ownedp;
2393 IS edgeNumberingIS_subdomain_new2old;
2394 ISCreateGeneral(PROTEUS_COMM_WORLD,edges_subdomain_owned.size(),&edgeNumbering_new2old[0],PETSC_COPY_VALUES,&edgeNumberingIS_subdomain_new2old);
2395 IS edgeNumberingIS_global_new2old;
2396 ISAllGather(edgeNumberingIS_subdomain_new2old,&edgeNumberingIS_global_new2old);
2397 const PetscInt *edgeNumbering_global_new2old;
2398 valarray<int> edgeNumbering_global_old2new(newMesh.
nEdges_global);
2399 ISGetIndices(edgeNumberingIS_global_new2old,&edgeNumbering_global_new2old);
2400 ierr =
enforceMemoryLimit(PROTEUS_COMM_WORLD, rank, max_rss_gb,
"Setting edgeNumering old2new/new2old");CHKERRABORT(PROTEUS_COMM_WORLD, ierr);
2403 edgeNumbering_global_old2new[edgeNumbering_global_new2old[edN]] = edN;
2405 ISRestoreIndices(edgeNumberingIS_global_new2old,&edgeNumbering_global_new2old);
2406 ISDestroy(&edgeNumberingIS_subdomain_new2old);
2407 ISDestroy(&edgeNumberingIS_global_new2old);
2408 ierr =
enforceMemoryLimit(PROTEUS_COMM_WORLD, rank, max_rss_gb,
"Done allocating edgeNumering old2new/new2old");CHKERRABORT(PROTEUS_COMM_WORLD, ierr);
2413 set<int> elements_overlap,nodes_overlap,elementBoundaries_overlap,edges_overlap;
2414 for (
int nN = 0; nN < nNodes_subdomain_new[rank]; nN++)
2417 for (
int offset = nodeStarOffsetsNew[nN];offset<nodeStarOffsetsNew[nN+1];offset++)
2419 int nN_point_global = nodeStarArrayNew[offset];
2420 bool offproc = nN_point_global < nodeOffsets_new[rank] || nN_point_global >= nodeOffsets_new[rank+1];
2422 nodes_overlap.insert(nN_point_global);
2425 for (
int eN_star_offset = nodeElementOffsets[nN];
2426 eN_star_offset < nodeElementOffsets[nN+1]; eN_star_offset++)
2428 int eN_star_old = nodeElementsArray[eN_star_offset];
2429 int eN_star_new = elementNumbering_global_old2new[eN_star_old];
2430 bool offproc = eN_star_new >= elementOffsets_new[rank+1] || eN_star_new < elementOffsets_new[rank];
2432 elements_overlap.insert(eN_star_new);
2435 for (
int ebN_star_offset = nodeElementBoundaryOffsets[nN];
2436 ebN_star_offset < nodeElementBoundaryOffsets[nN+1]; ebN_star_offset++)
2438 int ebN_star_old = nodeElementBoundariesArray[ebN_star_offset];
2439 int ebN_star_new = elementBoundaryNumbering_global_old2new[ebN_star_old];
2440 bool offproc = ebN_star_new >= elementBoundaryOffsets_new[rank+1] || ebN_star_new < elementBoundaryOffsets_new[rank];
2442 elementBoundaries_overlap.insert(ebN_star_new);
2445 for (
int edN_star_offset = nodeEdgeOffsets[nN];
2446 edN_star_offset < nodeEdgeOffsets[nN+1]; edN_star_offset++)
2448 int edN_star_old = nodeEdgesArray[edN_star_offset];
2449 int edN_star_new = edgeNumbering_global_old2new[edN_star_old];
2450 bool offproc = edN_star_new >= edgeOffsets_new[rank+1] || edN_star_new < edgeOffsets_new[rank];
2452 edges_overlap.insert(edN_star_new);
2455 elementNumbering_global_old2new.resize(0);
2457 assert(edges_overlap.size() + nEdges_subdomain_new[rank] == edgeNodesMap.size());
2461 int nN_subdomain = nNodes_subdomain_new[rank];
2462 map<int,int> nodes_overlap_global2subdomainMap;
2463 for (set<int>::iterator nN_globalp=nodes_overlap.begin();nN_globalp != nodes_overlap.end(); nN_globalp++,nN_subdomain++)
2464 nodes_overlap_global2subdomainMap[*nN_globalp] = nN_subdomain;
2483 PetscLogEventEnd(build_subdomains_edges_event,0,0,0,0);
2484 int build_subdomains_renumber_event;
2485 PetscLogEventRegister(
"Subd's renumber",0,&build_subdomains_renumber_event);
2486 PetscLogEventBegin(build_subdomains_renumber_event,0,0,0,0);
2491 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;
2492 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;
2510 map<int,int> nodeNumbering_global2subdomainMap;
2511 map<int,int> elementBoundaryNumbering_global2subdomainMap;
2512 map<int,int> edgeNumbering_global2subdomainMap;
2518 for (
int iv = 0; iv < nNodes_global; iv++)
2520 int nv;
double x,y,
z;
int nodeId(0);
2522 if (hasVertexMarkers > 0)
2523 vertexFile >> nodeId;
2525 assert(0 <= nv && nv < nNodes_global && vertexFile.good());
2526 int nN_global_new = nodeNumbering_global_old2new[nv];
2528 if (nN_global_new >= nodeOffsets_new[rank] && nN_global_new < nodeOffsets_new[rank+1])
2530 int nv_subdomain_new = nN_global_new - nodeOffsets_new[rank];
2531 nodeNumbering_subdomain2global[nv_subdomain_new] = nN_global_new;
2532 nodeNumbering_global2subdomainMap[nN_global_new] = nv_subdomain_new;
2536 if (hasVertexMarkers > 0)
2540 if (nodes_overlap.count(nN_global_new) == 1)
2542 int nv_subdomain_new = nodes_overlap_global2subdomainMap[nN_global_new];
2543 nodeNumbering_subdomain2global[nv_subdomain_new] = nN_global_new;
2544 nodeNumbering_global2subdomainMap[nN_global_new] = nv_subdomain_new;
2548 if (hasVertexMarkers > 0)
2554 ISRestoreIndices(nodeNumberingIS_global_old2new,&nodeNumbering_global_old2new);
2555 ISDestroy(&nodePartitioningIS_new);
2556 ISDestroy(&nodeNumberingIS_subdomain_old2new);
2557 ISDestroy(&nodeNumberingIS_global_old2new);
2559 ierr =
enforceMemoryLimit(PROTEUS_COMM_WORLD, rank, max_rss_gb,
"Done reading vertices");CHKERRABORT(PROTEUS_COMM_WORLD, ierr);
2568 for (
int eN = 0; eN < nElements_subdomain_new[rank]; eN++)
2570 int eN_global_new = elementOffsets_new[rank] + eN;
2571 int eN_global_old = elementNumbering_global_new2old[eN_global_new];
2572 elementNumbering_subdomain2global[eN] = eN_global_new;
2576 int nN_global_new = elementNodesArrayMap[eN_global_old][nN];
2577 int nN_subdomain = nodeNumbering_global2subdomainMap[nN_global_new];
2584 set<int>::iterator eN_p = elements_overlap.begin();
2585 for (
int eN = nElements_subdomain_new[rank]; eN < nElements_subdomain_new[rank] + int(elements_overlap.size()); eN++,eN_p++)
2587 int eN_global_new = *eN_p;
2588 int eN_global_old = elementNumbering_global_new2old[eN_global_new];
2589 elementNumbering_subdomain2global[eN] = eN_global_new;
2593 int nN_global_new = elementNodesArrayMap[eN_global_old][nN];
2594 int nN_subdomain = nodeNumbering_global2subdomainMap[nN_global_new];
2598 ISRestoreIndices(elementNumberingIS_global_new2old,&elementNumbering_global_new2old);
2599 ISDestroy(&elementNumberingIS_subdomain_new2old);
2600 ISDestroy(&elementNumberingIS_global_new2old);
2606 for (
int ebN=0; ebN < nElementBoundaries_subdomain_new[rank]; ebN++)
2608 int ebN_global = ebN + elementBoundaryOffsets_new[rank];
2609 elementBoundaryNumbering_subdomain2global[ebN]=ebN_global;
2610 elementBoundaryNumbering_global2subdomainMap[ebN_global] = ebN;
2615 set<int>::iterator ebN_p = elementBoundaries_overlap.begin();
2616 for(
int ebN=nElementBoundaries_subdomain_new[rank];ebN < nElementBoundaries_subdomain_new[rank] + int(elementBoundaries_overlap.size()); ebN++,ebN_p++)
2618 int ebN_global = *ebN_p;
2619 elementBoundaryNumbering_subdomain2global[ebN] = ebN_global;
2620 elementBoundaryNumbering_global2subdomainMap[ebN_global] = ebN;
2629 for (
int eN=0;eN<nElements_subdomain_new[rank];eN++)
2631 int eN_global = eN+elementOffsets_new[rank];
2635 elementBoundaryNumbering_global2subdomainMap[elementBoundaryNumbering_global_old2new[elementBoundariesMap[eN_global][ebN]]];
2641 eN_p = elements_overlap.begin();
2642 for (
int eN = nElements_subdomain_new[rank]; eN < nElements_subdomain_new[rank] + int(elements_overlap.size()); eN++,eN_p++)
2644 int eN_global_new = *eN_p;
2648 elementBoundaryNumbering_global2subdomainMap[elementBoundaryNumbering_global_old2new[elementBoundariesMap[eN_global_new][ebN]]];
2656 for (
int edN=0; edN < nEdges_subdomain_new[rank]; edN++)
2658 int edN_global = edN + edgeOffsets_new[rank];
2659 edgeNumbering_subdomain2global[edN]=edN_global;
2660 edgeNumbering_global2subdomainMap[edN_global] = edN;
2665 set<int>::iterator edN_p = edges_overlap.begin();
2666 for(
int edN=nEdges_subdomain_new[rank];edN < nEdges_subdomain_new[rank] + int(edges_overlap.size()); edN++,edN_p++)
2668 int edN_global = *edN_p;
2669 edgeNumbering_subdomain2global[edN] = edN_global;
2670 edgeNumbering_global2subdomainMap[edN_global] = edN;
2678 for (map<
int,pair<int,int> >::iterator edgep=edgeNodesMap.begin();
2679 edgep!=edgeNodesMap.end();
2682 int edN_global_old = edgep->first;
2683 int edN_global_new = edgeNumbering_global_old2new[edN_global_old];
2684 assert(edgeNumbering_global2subdomainMap.find(edN_global_new) != edgeNumbering_global2subdomainMap.end());
2685 int edN_subdomain = edgeNumbering_global2subdomainMap[edN_global_new];
2689 edgeNumbering_global_old2new.resize(0);
2710 using namespace std;
2715 int> elementBoundaryIds;
2722 register int nodes[3];
2727 if(elementBoundaryElements.find(ebt) != elementBoundaryElements.end())
2729 elementBoundaryElements[ebt].right=eN;
2730 elementBoundaryElements[ebt].right_ebN_element=ebN;
2731 assert(elementBoundaryIds[ebt] == ebN_global);
2735 elementBoundaryElements.insert(elementBoundaryElements.end(),make_pair(ebt,
ElementNeighbors(eN,ebN)));
2736 elementBoundaryIds.insert(elementBoundaryIds.end(),make_pair(ebt,ebN_global));
2746 set<int> interiorElementBoundaries,exteriorElementBoundaries;
2757 eb != elementBoundaryElements.end();
2760 int ebN = elementBoundaryIds[eb->first];
2770 if(eb->second.right != -1)
2772 interiorElementBoundaries.insert(ebN);
2776 exteriorElementBoundaries.insert(ebN);
2778 if (eb->second.right != -1)
2788 for (set<int>::iterator ebN=interiorElementBoundaries.begin();ebN != interiorElementBoundaries.end(); ebN++,ebNI++)
2790 for (set<int>::iterator ebN=exteriorElementBoundaries.begin();ebN != exteriorElementBoundaries.end(); ebN++,ebNE++)
2792 set<NodeTuple<2> > edges;
2817 for (set<int>::iterator nN_star=nodeStar[nN].begin();nN_star!=nodeStar[nN].end();nN_star++,offset++)
2837 for (set<int>::iterator eN_star = nodeElementsStar[nN].begin(); eN_star != nodeElementsStar[nN].end();
2876 if (hasElementBoundaryMarkers)
2879 for (map<int,int>::iterator ebmp = elementBoundaryMaterialTypesMap.begin(); ebmp != elementBoundaryMaterialTypesMap.end();ebmp++)
2881 int ebN_global_new = elementBoundaryNumbering_global_old2new[ebmp->first];
2882 assert(elementBoundaryNumbering_global2subdomainMap.find(ebN_global_new) != elementBoundaryNumbering_global2subdomainMap.end());
2883 int ebN_subdomain = elementBoundaryNumbering_global2subdomainMap[ebN_global_new];
2886 if (!hasVertexMarkers)
2897 elementBoundaryNumbering_global_old2new.resize(0);
2898 ierr =
enforceMemoryLimit(PROTEUS_COMM_WORLD, rank, max_rss_gb,
"Done with material types");CHKERRABORT(PROTEUS_COMM_WORLD, ierr);
2899 PetscLogEventEnd(build_subdomains_renumber_event,0,0,0,0);
2900 int build_subdomains_cleanup_event;
2901 PetscLogEventRegister(
"Cleanup",0,&build_subdomains_cleanup_event);
2902 PetscLogEventBegin(build_subdomains_cleanup_event,0,0,0,0);
2916 for (
int sdN = 0; sdN < size+1; sdN++)
2947 H5Sclose(filespace);
2952 PetscLogEventEnd(build_subdomains_cleanup_event,0,0,0,0);
2954 PetscLogView(PETSC_VIEWER_STDOUT_WORLD);
2955 ierr =
enforceMemoryLimit(PROTEUS_COMM_WORLD, rank, max_rss_gb,
"Done with partitioning!");CHKERRABORT(PROTEUS_COMM_WORLD, ierr);
2960 using namespace std;
2961 PetscErrorCode ierr;
2962 PetscMPIInt size,rank;
2964 ierr = MPI_Comm_size(PROTEUS_COMM_WORLD,&size);CHKERRABORT(PROTEUS_COMM_WORLD, ierr);
2965 ierr = MPI_Comm_rank(PROTEUS_COMM_WORLD,&rank);CHKERRABORT(PROTEUS_COMM_WORLD, ierr);
2966 PetscLogStage partitioning_stage;
2967 PetscLogStageRegister(
"Mesh Partition",&partitioning_stage);
2968 PetscLogStagePush(partitioning_stage);
2995 bool failed =
false;
2996 const int simplexDim = 3;
2997 const int vertexDim = 3;
2999 std::string vertexFileName = std::string(filebase) +
".node" ;
3000 std::string elementFileName = std::string(filebase) +
".ele" ;
3001 std::string elementBoundaryFileName = std::string(filebase) +
".edge" ;
3002 std::string edgeFileName = std::string(filebase) +
".edge" ;
3013 int read_elements_event;
3014 PetscLogEventRegister(
"Read eles",0,&read_elements_event);
3015 PetscLogEventBegin(read_elements_event,0,0,0,0);
3016 std::ifstream vertexFile(vertexFileName.c_str());
3017 if (!vertexFile.good())
3019 std::cerr<<
"cannot open triangle node file "
3020 <<vertexFileName<<std::endl;
3024 int hasVertexMarkers(0),hasVertexAttributes(0),nSpace(2),nNodes_global;
3026 vertexFile >>
eatcomments >> nNodes_global >> nSpace >> hasVertexAttributes >> hasVertexMarkers >>
eatline ;
3027 assert(nNodes_global > 0);
3028 assert(nSpace == 2);
3033 if (hasVertexAttributes > 0)
3035 std::cerr<<
"WARNING triangle nodes hasAttributes= "<<hasVertexAttributes
3036 <<
" > 0 will treat first value as integer id for boundary!!"<<std::endl;
3037 hasVertexMarkers = 1;
3043 valarray<int> nodeOffsets_old(size+1);
3044 nodeOffsets_old[0] = 0;
3045 for (
int sdN=0; sdN < size; sdN++)
3047 nodeOffsets_old[sdN+1] = nodeOffsets_old[sdN] +
3048 int(nNodes_global)/size + (int(nNodes_global)%size > sdN);
3050 int nNodes_subdomain_old = nodeOffsets_old[rank+1] - nodeOffsets_old[rank];
3057 std::ifstream elementFile(elementFileName.c_str());
3058 if (!elementFile.good())
3060 std::cerr<<
"cannot open triangle element file "
3061 <<elementFileName<<std::endl;
3066 int nNodesPerSimplex(simplexDim),hasElementMarkers = 0,nElements_global;
3067 elementFile >>
eatcomments >> nElements_global >> nNodesPerSimplex >> hasElementMarkers >>
eatline;
3068 assert(nElements_global > 0);
3069 assert(nNodesPerSimplex == simplexDim);
3071 vector<int> element_nodes_old(3);
3072 vector<set<int> > nodeStar(nNodes_subdomain_old);
3073 map<int,vector<int> > elements_old;
3074 for (
int ie = 0; ie < nElements_global; ie++)
3079 assert(0 <= ne && ne < nElements_global && elementFile.good());
3080 for (
int iv = 0; iv < simplexDim; iv++)
3084 assert(0 <= nv && nv < nNodes_global);
3085 element_nodes_old[iv] = nv;
3088 for (
int iv = 0; iv < simplexDim; iv++)
3091 int nN_star = element_nodes_old[iv];
3092 bool inSubdomain=
false;
3093 if (nN_star >= nodeOffsets_old[rank] && nN_star < nodeOffsets_old[rank+1])
3097 for (
int jv = 0; jv < simplexDim; jv++)
3101 int nN_star_subdomain = nN_star-nodeOffsets_old[rank];
3102 nodeStar[nN_star_subdomain].insert(element_nodes_old[jv]);
3107 elements_old[ie] = element_nodes_old;
3111 elementFile.close();
3112 PetscLogEventEnd(read_elements_event,0,0,0,0);
3113 int repartition_nodes_event;
3114 PetscLogEventRegister(
"Repart nodes",0,&repartition_nodes_event);
3115 PetscLogEventBegin(repartition_nodes_event,0,0,0,0);
3118 valarray<int> nodeStarOffsets(nNodes_subdomain_old+1);
3119 nodeStarOffsets[0] = 0;
3120 for (
int nN=1;nN<nNodes_subdomain_old+1;nN++)
3121 nodeStarOffsets[nN] = nodeStarOffsets[nN-1] + nodeStar[nN-1].size();
3122 valarray<int> nodeStarArray(nodeStarOffsets[nNodes_subdomain_old]);
3123 for (
int nN=0,offset=0;nN<nNodes_subdomain_old;nN++)
3124 for (set<int>::iterator nN_star=nodeStar[nN].begin();nN_star!=nodeStar[nN].end();nN_star++,offset++)
3125 nodeStarArray[offset] = *nN_star;
3127 int max_nNodeNeighbors_node=0;
3128 for (
int nN=0;nN<nNodes_subdomain_old;nN++)
3129 max_nNodeNeighbors_node=
max(max_nNodeNeighbors_node,nodeStarOffsets[nN+1]-nodeStarOffsets[nN]);
3131 PetscBool isInitialized;
3132 PetscInitialized(&isInitialized);
3133 PetscInt *nodeNeighborsOffsets_subdomain,*nodeNeighbors_subdomain,*weights_subdomain,*vertex_weights_subdomain;
3134 PetscReal *partition_weights;
3135 PetscMalloc(
sizeof(PetscInt)*(nNodes_subdomain_old+1),&nodeNeighborsOffsets_subdomain);
3136 PetscMalloc(
sizeof(PetscInt)*(nNodes_subdomain_old*max_nNodeNeighbors_node),&nodeNeighbors_subdomain);
3137 PetscMalloc(
sizeof(PetscInt)*(nNodes_subdomain_old*max_nNodeNeighbors_node),&weights_subdomain);
3138 PetscMalloc(
sizeof(PetscInt)*(nNodes_subdomain_old),&vertex_weights_subdomain);
3139 PetscMalloc(
sizeof(PetscReal)*(size),&partition_weights);
3140 for (
int sd=0;sd<size;sd++)
3141 partition_weights[sd] = 1.0/
double(size);
3142 nodeNeighborsOffsets_subdomain[0] = 0;
3144 for (
int nN = 0,offset=0; nN < nNodes_subdomain_old; nN++)
3146 for (
int offset_subdomain = nodeStarOffsets[nN];
3147 offset_subdomain < nodeStarOffsets[nN+1];
3150 nodeNeighbors_subdomain[offset++] = nodeStarArray[offset_subdomain];
3152 nodeNeighborsOffsets_subdomain[nN+1]=offset;
3153 sort(&nodeNeighbors_subdomain[nodeNeighborsOffsets_subdomain[nN]],&nodeNeighbors_subdomain[nodeNeighborsOffsets_subdomain[nN+1]]);
3155 int weight= (nodeNeighborsOffsets_subdomain[nN+1] - nodeNeighborsOffsets_subdomain[nN]);
3156 vertex_weights_subdomain[nN] = weight;
3157 for (
int k=nodeNeighborsOffsets_subdomain[nN];k<nodeNeighborsOffsets_subdomain[nN+1];k++)
3158 weights_subdomain[k] = weight;
3164 int nNodes_subdomain_max=0;
3165 MPI_Allreduce(&nNodes_subdomain_old,
3166 &nNodes_subdomain_max,
3170 PROTEUS_COMM_WORLD);
3172 std::cout<<
"Max nNodes_subdomain "<<nNodes_subdomain_max<<
" nNodes_global "<<nNodes_global<<std::endl;
3173 ierr = MatCreateMPIAdj(PROTEUS_COMM_WORLD,
3174 nNodes_subdomain_old,
3176 nodeNeighborsOffsets_subdomain,
3177 nodeNeighbors_subdomain,
3179 &petscAdjacency);CHKERRABORT(PROTEUS_COMM_WORLD, ierr);
3181 const double max_rss_gb(3.75);
3182 ierr =
enforceMemoryLimit(PROTEUS_COMM_WORLD, rank, max_rss_gb,
"Done allocating MPIAdj");CHKERRABORT(PROTEUS_COMM_WORLD, ierr);
3183 MatPartitioning petscPartition;
3184 ierr = MatPartitioningCreate(PROTEUS_COMM_WORLD,&petscPartition);CHKERRABORT(PROTEUS_COMM_WORLD, ierr);
3185 ierr = MatPartitioningSetAdjacency(petscPartition,petscAdjacency);CHKERRABORT(PROTEUS_COMM_WORLD, ierr);
3186 ierr = MatPartitioningSetFromOptions(petscPartition);CHKERRABORT(PROTEUS_COMM_WORLD, ierr);
3187 ierr = MatPartitioningSetVertexWeights(petscPartition,vertex_weights_subdomain);CHKERRABORT(PROTEUS_COMM_WORLD, ierr);
3188 ierr = MatPartitioningSetPartitionWeights(petscPartition,partition_weights);CHKERRABORT(PROTEUS_COMM_WORLD, ierr);
3190 IS nodePartitioningIS_new;
3191 ierr = MatPartitioningApply(petscPartition,&nodePartitioningIS_new);CHKERRABORT(PROTEUS_COMM_WORLD, ierr);
3192 ierr = MatPartitioningDestroy(&petscPartition);CHKERRABORT(PROTEUS_COMM_WORLD, ierr);
3193 ierr =
enforceMemoryLimit(PROTEUS_COMM_WORLD, rank, max_rss_gb,
"Done applying partition");CHKERRABORT(PROTEUS_COMM_WORLD, ierr);
3196 valarray<int> nNodes_subdomain_new(size);
3197 ISPartitioningCount(nodePartitioningIS_new,size,&nNodes_subdomain_new[0]);
3200 valarray<int> nodeOffsets_new(size+1);
3201 nodeOffsets_new[0] = 0;
3202 for (
int sdN = 0; sdN < size; sdN++)
3204 nodeOffsets_new[sdN+1] = nodeOffsets_new[sdN] + nNodes_subdomain_new[sdN];
3208 IS nodeNumberingIS_subdomain_old2new;
3209 ISPartitioningToNumbering(nodePartitioningIS_new,&nodeNumberingIS_subdomain_old2new);
3216 MPI_Info info = MPI_INFO_NULL;
3217 hid_t plist_id = H5Pcreate(H5P_FILE_ACCESS);
3218 H5Pset_fapl_mpio(plist_id, PROTEUS_COMM_WORLD, info);
3223 const char* H5FILE_NAME(
"mappings.h5");
3224 hid_t file_id = H5Fcreate(H5FILE_NAME, H5F_ACC_TRUNC, H5P_DEFAULT, plist_id);
3232 dimsf[0] = nNodes_global;
3234 hid_t filespace = H5Screate_simple(
RANK, dimsf, NULL);
3239 hid_t dset_id = H5Dcreate(file_id,
"nodeNumbering_old2new", H5T_NATIVE_INT, filespace,
3240 H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
3241 H5Sclose(filespace);
3249 count[0] = nNodes_subdomain_old;
3250 offset[0] = nodeOffsets_old[rank];
3251 hid_t memspace = H5Screate_simple(
RANK, count, NULL);
3256 filespace = H5Dget_space(dset_id);
3257 H5Sselect_hyperslab(filespace, H5S_SELECT_SET, offset, NULL, count, NULL);
3266 const PetscInt* data;
3267 ISGetIndices(nodeNumberingIS_subdomain_old2new, &data);
3272 plist_id = H5Pcreate(H5P_DATASET_XFER);
3273 H5Pset_dxpl_mpio(plist_id, H5FD_MPIO_COLLECTIVE);
3275 herr_t status = H5Dwrite(dset_id, H5T_NATIVE_INT, memspace, filespace,
3278 ISRestoreIndices(nodeNumberingIS_subdomain_old2new, &data);
3289 IS nodeNumberingIS_global_old2new;
3290 ISAllGather(nodeNumberingIS_subdomain_old2new,&nodeNumberingIS_global_old2new);
3291 const PetscInt * nodeNumbering_global_old2new;
3292 ISGetIndices(nodeNumberingIS_global_old2new,&nodeNumbering_global_old2new);
3300 int dset_data[nNodes_global];
3306 dataset_id = H5Dopen2(file_id,
"/nodeNumbering_old2new", H5P_DEFAULT);
3308 status = H5Dread(dataset_id, H5T_NATIVE_INT, H5S_ALL, H5S_ALL, H5P_DEFAULT,
3312 status = H5Dclose(dataset_id);
3314 for (
int i=0;i<nNodes_global;i++)
3315 assert(nodeNumbering_global_old2new[i] == dset_data[i]);
3316 std::cout<<
"==================out of core old2new is correct!===================="<<std::endl;
3328 PetscLogEventEnd(repartition_nodes_event,0,0,0,0);
3329 int receive_element_mask_event;
3330 PetscLogEventRegister(
"Recv. ele mask",0,&receive_element_mask_event);
3331 PetscLogEventBegin(receive_element_mask_event,0,0,0,0);
3336 PetscLogEventEnd(receive_element_mask_event,0,0,0,0);
3337 ierr =
enforceMemoryLimit(PROTEUS_COMM_WORLD, rank, max_rss_gb,
"Done with masks");CHKERRABORT(PROTEUS_COMM_WORLD, ierr);
3338 int build_subdomains_reread_elements_event;
3339 PetscLogEventRegister(
"Reread eles",0,&build_subdomains_reread_elements_event);
3340 PetscLogEventBegin(build_subdomains_reread_elements_event,0,0,0,0);
3347 std::ifstream elementFile2(elementFileName.c_str());
3349 if (!elementFile2.good())
3351 std::cerr<<
"cannot open triangle elements file"
3352 <<elementFileName<<std::endl;
3356 elementFile2 >>
eatcomments >> nElements_global >> nNodesPerSimplex >> hasElementMarkers >>
eatline;
3357 assert(nElements_global > 0);
3358 assert(nNodesPerSimplex == simplexDim);
3359 set<int> elements_subdomain_owned;
3360 vector<int> element_nodes_new(3);
3361 int element_nodes_new_array[3];
3362 vector<set<int> > nodeElementsStar(nNodes_subdomain_new[rank]);
3363 vector<set<int> > nodeStarNew(nNodes_subdomain_new[rank]);
3364 map<int,vector<int> > elementNodesArrayMap;
3365 map<int,long int> elementMaterialTypesMap;
3367 map<NodeTuple<2>,set<pair<int,int> > > edgeElementsMap;
3369 for (
int ie = 0; ie < nElements_global; ie++)
3371 int ne, nv, elementId(0);
3372 long double elementId_double;
3375 assert(0 <= ne && ne < nElements_global && elementFile.good());
3376 for (
int iv = 0; iv < simplexDim; iv++)
3378 elementFile2 >> nv ;
3380 assert(0 <= nv && nv < nNodes_global);
3381 element_nodes_old[iv] = nv;
3382 element_nodes_new[iv] = nodeNumbering_global_old2new[nv];
3383 element_nodes_new_array[iv] = element_nodes_new[iv];
3386 for (
int iv = 0; iv < simplexDim; iv++)
3388 int nN_star_new = element_nodes_new[iv];
3389 bool inSubdomain=
false;
3390 if (nN_star_new >= nodeOffsets_new[rank] && nN_star_new < nodeOffsets_new[rank+1])
3394 for (
int ebN=0;ebN < simplexDim ; ebN++)
3396 int nodes[simplexDim-1] = { element_nodes_new[(ebN+1) % simplexDim],
3397 element_nodes_new[(ebN+2) % simplexDim]};
3398 NodeTuple<simplexDim-1> nodeTuple(nodes);
3399 if(elementBoundaryElementsMap.find(nodeTuple) != elementBoundaryElementsMap.end())
3401 if (elementBoundaryElementsMap[nodeTuple].right == -1 && ne != elementBoundaryElementsMap[nodeTuple].left)
3403 elementBoundaryElementsMap[nodeTuple].right=ne;
3404 elementBoundaryElementsMap[nodeTuple].right_ebN_element=ebN;
3413 for (
int nNL=0,edN=0;nNL < simplexDim ; nNL++)
3414 for(
int nNR=nNL+1;nNR < simplexDim;nNR++,edN++)
3416 int nodes[2] = { element_nodes_new[nNL],
3417 element_nodes_new[nNR]};
3419 edgeElementsMap[nodeTuple].insert(pair<int,int>(ne,edN));
3422 int nN_star_new_subdomain = nN_star_new - nodeOffsets_new[rank];
3423 nodeElementsStar[nN_star_new_subdomain].insert(ne);
3424 for (
int jv = 0; jv < simplexDim; jv++)
3428 int nN_point_new = element_nodes_new[jv];
3429 nodeStarNew[nN_star_new_subdomain].insert(nN_point_new);
3435 elementNodesArrayMap[ne] = element_nodes_new;
3438 if (elementNodesArrayMap.find(ne) != elementNodesArrayMap.end())
3440 if (nodeTuple.
nodes[1] >= nodeOffsets_new[rank] && nodeTuple.
nodes[1] < nodeOffsets_new[rank+1])
3441 elements_subdomain_owned.insert(ne);
3442 if (hasElementMarkers > 0)
3444 elementFile2 >> elementId_double;
3445 elementId =
static_cast<long int>(elementId_double);
3446 elementMaterialTypesMap[ne] = elementId;
3451 elementFile2.close();
3452 int nElements_owned_subdomain(elements_subdomain_owned.size()),
3453 nElements_owned_new=0;
3454 MPI_Allreduce(&nElements_owned_subdomain,&nElements_owned_new,1,MPI_INT,MPI_SUM,PROTEUS_COMM_WORLD);
3455 assert(nElements_owned_new == nElements_global);
3456 PetscLogEventEnd(build_subdomains_reread_elements_event,0,0,0,0);
3457 int build_subdomains_send_marked_elements_event;
3458 PetscLogEventRegister(
"Mark/send eles",0,&build_subdomains_send_marked_elements_event);
3459 PetscLogEventBegin(build_subdomains_send_marked_elements_event,0,0,0,0);
3460 ierr =
enforceMemoryLimit(PROTEUS_COMM_WORLD, rank, max_rss_gb,
"Done marking elements");CHKERRABORT(PROTEUS_COMM_WORLD, ierr);
3465 valarray<int> nodeElementOffsets(nNodes_subdomain_new[rank]+1);
3466 nodeElementOffsets[0] = 0;
3467 for (
int nN = 0; nN < nNodes_subdomain_new[rank]; nN++)
3468 nodeElementOffsets[nN+1] = nodeElementOffsets[nN]+nodeElementsStar[nN].size();
3469 valarray<int> nodeElementsArray(nodeElementOffsets[nNodes_subdomain_new[rank]]);
3470 for (
int nN=0,offset=0; nN < nNodes_subdomain_new[rank]; nN++)
3472 for (set<int>::iterator eN_star = nodeElementsStar[nN].begin(); eN_star != nodeElementsStar[nN].end();
3475 nodeElementsArray[offset] = *eN_star;
3479 valarray<int> nodeStarOffsetsNew(nNodes_subdomain_new[rank]+1);
3480 nodeStarOffsetsNew[0] = 0;
3481 for (
int nN=1;nN<nNodes_subdomain_new[rank]+1;nN++)
3482 nodeStarOffsetsNew[nN] = nodeStarOffsetsNew[nN-1] + nodeStarNew[nN-1].size();
3483 valarray<int> nodeStarArrayNew(nodeStarOffsetsNew[nNodes_subdomain_new[rank]]);
3484 for (
int nN=0,offset=0;nN<nNodes_subdomain_new[rank];nN++)
3486 for (set<int>::iterator nN_star=nodeStarNew[nN].begin();nN_star!=nodeStarNew[nN].end();nN_star++,offset++)
3488 nodeStarArrayNew[offset] = *nN_star;
3491 PetscLogEventEnd(build_subdomains_send_marked_elements_event,0,0,0,0);
3492 int build_subdomains_global_numbering_elements_event;
3493 PetscLogEventRegister(
"Global ele nmbr",0,&build_subdomains_global_numbering_elements_event);
3494 PetscLogEventBegin(build_subdomains_global_numbering_elements_event,0,0,0,0);
3498 valarray<int> nElements_subdomain_new(size),
3499 elementOffsets_new(size+1);
3500 for (
int sdN = 0; sdN < size; sdN++)
3503 nElements_subdomain_new[sdN] = int(elements_subdomain_owned.size());
3505 nElements_subdomain_new[sdN] = 0;
3507 valarray<int> nElements_subdomain_new_send = nElements_subdomain_new;
3508 MPI_Allreduce(&nElements_subdomain_new_send[0],&nElements_subdomain_new[0],size,MPI_INT,MPI_SUM,PROTEUS_COMM_WORLD);
3510 elementOffsets_new[0] = 0;
3511 for (
int sdN = 0; sdN < size; sdN++)
3512 elementOffsets_new[sdN+1] = elementOffsets_new[sdN] + nElements_subdomain_new[sdN];
3514 valarray<int> elementNumbering_subdomain_new2old(elements_subdomain_owned.size());
3515 set<int>::iterator eN_ownedp = elements_subdomain_owned.begin();
3516 for (
int eN = 0; eN < int(elements_subdomain_owned.size()); eN++,eN_ownedp++)
3518 elementNumbering_subdomain_new2old[eN] = *eN_ownedp;
3521 IS elementNumberingIS_subdomain_new2old;
3522 ISCreateGeneral(PROTEUS_COMM_WORLD,elements_subdomain_owned.size(),&elementNumbering_subdomain_new2old[0],PETSC_COPY_VALUES,
3523 &elementNumberingIS_subdomain_new2old);
3524 IS elementNumberingIS_global_new2old;
3525 ISAllGather(elementNumberingIS_subdomain_new2old,&elementNumberingIS_global_new2old);
3527 const PetscInt *elementNumbering_global_new2old;
3528 ISGetIndices(elementNumberingIS_global_new2old,&elementNumbering_global_new2old);
3530 valarray<int> elementNumbering_global_old2new(nElements_global);
3531 for (
int eN = 0; eN < nElements_global; eN++)
3533 elementNumbering_global_old2new[elementNumbering_global_new2old[eN]] = eN;
3535 PetscLogEventEnd(build_subdomains_global_numbering_elements_event,0,0,0,0);
3536 ierr =
enforceMemoryLimit(PROTEUS_COMM_WORLD, rank, max_rss_gb,
"Done allocating element numbering new2old/old2new");CHKERRABORT(PROTEUS_COMM_WORLD, ierr);
3537 int build_subdomains_faces_event;
3538 PetscLogEventRegister(
"Subd faces",0,&build_subdomains_faces_event);
3539 PetscLogEventBegin(build_subdomains_faces_event,0,0,0,0);
3551 std::ifstream elementBoundaryFile(elementBoundaryFileName.c_str());
3553 if (!elementBoundaryFile.good())
3555 std::cerr<<
"cannot open triangle edge file "
3556 <<elementBoundaryFileName<<std::endl;
3561 bool hasElementBoundaryMarkers =
false;
3562 int nElementBoundaries_global;
3563 int ihasElementBoundaryMarkers(0);
3565 elementBoundaryFile >>
eatcomments >> nElementBoundaries_global >>ihasElementBoundaryMarkers >>
eatline ;
3566 assert(nElementBoundaries_global > 0);
3567 if (ihasElementBoundaryMarkers > 0)
3569 hasElementBoundaryMarkers =
true;
3573 set<int> elementBoundaries_subdomain_owned;
3574 vector<set<int> > nodeElementBoundariesStar(nNodes_subdomain_new[rank]);
3575 map<int,int> elementBoundaryMaterialTypesMap;
3576 map<int,vector<int> > elementBoundariesMap;
3577 set<int> supportedElementBoundaries;
3578 for (
int ieb = 0; ieb < nElementBoundaries_global; ieb++)
3580 int neb,nn0,nn1;
int ebId(0);
3581 elementBoundaryFile >>
eatcomments >> neb >> nn0 >> nn1;
3582 if (ihasElementBoundaryMarkers > 0)
3583 elementBoundaryFile >> ebId;
3587 assert(0 <= neb && neb < nElementBoundaries_global && elementBoundaryFile.good());
3590 int nn0_new = nodeNumbering_global_old2new[nn0];
3591 if (nn0_new >= nodeOffsets_new[rank] && nn0_new < nodeOffsets_new[rank+1])
3593 nodeElementBoundariesStar[nn0_new-nodeOffsets_new[rank]].insert(neb);
3594 supportedElementBoundaries.insert(neb);
3596 int nn1_new = nodeNumbering_global_old2new[nn1];
3597 if (nn1_new >= nodeOffsets_new[rank] && nn1_new < nodeOffsets_new[rank+1])
3599 nodeElementBoundariesStar[nn1_new-nodeOffsets_new[rank]].insert(neb);
3600 supportedElementBoundaries.insert(neb);
3602 int nodes[2] = {nn0_new,nn1_new};
3604 elementBoundaryFile >>
eatline;
3605 if (elementBoundaryElementsMap.find(nodeTuple) != elementBoundaryElementsMap.end())
3607 if (nodeTuple.
nodes[1] >= nodeOffsets_new[rank] && nodeTuple.
nodes[1] < nodeOffsets_new[rank+1])
3608 elementBoundaries_subdomain_owned.insert(neb);
3609 if (ihasElementBoundaryMarkers > 0)
3610 elementBoundaryMaterialTypesMap[neb]=ebId;
3611 int eN_left = elementNumbering_global_old2new[elementBoundaryElementsMap[nodeTuple].left];
3612 if (elementBoundariesMap.find(eN_left) != elementBoundariesMap.end())
3614 elementBoundariesMap[eN_left][elementBoundaryElementsMap[nodeTuple].left_ebN_element] = neb;
3619 vector<int> elementBoundaries_element(3,-1);
3620 elementBoundariesMap[eN_left] = elementBoundaries_element;
3622 elementBoundariesMap[eN_left][elementBoundaryElementsMap[nodeTuple].left_ebN_element] = neb;
3624 if (elementBoundaryElementsMap[nodeTuple].right >= 0)
3626 int eN_right = elementNumbering_global_old2new[elementBoundaryElementsMap[nodeTuple].right];
3627 if (elementBoundariesMap.find(eN_right) != elementBoundariesMap.end())
3629 elementBoundariesMap[eN_right][elementBoundaryElementsMap[nodeTuple].right_ebN_element] = neb;
3634 vector<int> elementBoundaries_element(3,-1);
3635 elementBoundariesMap[eN_right] = elementBoundaries_element;
3637 elementBoundariesMap[eN_right][elementBoundaryElementsMap[nodeTuple].right_ebN_element] = neb;
3643 elementBoundaryFile.close();
3644 ierr =
enforceMemoryLimit(PROTEUS_COMM_WORLD, rank, max_rss_gb,
"Done reading element boundaries");CHKERRABORT(PROTEUS_COMM_WORLD, ierr);
3645 int nElementBoundaries_owned_subdomain=elementBoundaries_subdomain_owned.size(),
3646 nElementBoundaries_owned_new=0;
3647 MPI_Allreduce(&nElementBoundaries_owned_subdomain,&nElementBoundaries_owned_new,1,MPI_INT,MPI_SUM,PROTEUS_COMM_WORLD);
3648 assert(nElementBoundaries_owned_new == nElementBoundaries_global);
3651 for (map<
int,vector<int> >::iterator elementBoundariesp=elementBoundariesMap.begin();
3652 elementBoundariesp!=elementBoundariesMap.end();
3653 elementBoundariesp++)
3656 for (
int iv=0;iv<simplexDim;iv++)
3659 int nN_global = elementNodesArrayMap[elementNumbering_global_new2old[elementBoundariesp->first]][iv];
3660 if (nN_global >= nodeOffsets_new[rank] && nN_global < nodeOffsets_new[rank+1])
3663 for(
int eb=0;eb<simplexDim;eb++)
3665 nodeElementBoundariesStar[nN_global-nodeOffsets_new[rank]].insert(elementBoundariesp->second[eb]);
3671 valarray<int> nodeElementBoundaryOffsets(nNodes_subdomain_new[rank]+1);
3672 nodeElementBoundaryOffsets[0] = 0;
3673 for (
int nN = 0; nN < nNodes_subdomain_new[rank]; nN++)
3674 nodeElementBoundaryOffsets[nN+1] = nodeElementBoundaryOffsets[nN]+nodeElementBoundariesStar[nN].size();
3675 valarray<int> nodeElementBoundariesArray(nodeElementBoundaryOffsets[nNodes_subdomain_new[rank]]);
3676 for (
int nN=0,offset=0; nN < nNodes_subdomain_new[rank]; nN++)
3678 for (set<int>::iterator ebN_star = nodeElementBoundariesStar[nN].begin(); ebN_star != nodeElementBoundariesStar[nN].end();
3679 ebN_star++,offset++)
3681 nodeElementBoundariesArray[offset] = *ebN_star;
3686 valarray<int> nElementBoundaries_subdomain_new(size),
3687 elementBoundaryOffsets_new(size+1);
3688 for (
int sdN=0;sdN<size;sdN++)
3690 nElementBoundaries_subdomain_new[sdN] = elementBoundaries_subdomain_owned.size();
3692 nElementBoundaries_subdomain_new[sdN] = 0;
3693 valarray<int> nElementBoundaries_subdomain_new_send=nElementBoundaries_subdomain_new;
3694 MPI_Allreduce(&nElementBoundaries_subdomain_new_send[0],&nElementBoundaries_subdomain_new[0],size,MPI_INT,MPI_SUM,PROTEUS_COMM_WORLD);
3695 elementBoundaryOffsets_new[0] = 0;
3696 for (
int sdN=0;sdN<size;sdN++)
3697 elementBoundaryOffsets_new[sdN+1] = elementBoundaryOffsets_new[sdN]+nElementBoundaries_subdomain_new[sdN];
3703 valarray<int> elementBoundaryNumbering_new2old(elementBoundaries_subdomain_owned.size());
3704 set<int>::iterator ebN_ownedp=elementBoundaries_subdomain_owned.begin();
3705 for (
int ebN=0;ebN<int(elementBoundaries_subdomain_owned.size());ebN++)
3707 elementBoundaryNumbering_new2old[ebN] = *ebN_ownedp++;
3709 IS elementBoundaryNumberingIS_subdomain_new2old;
3710 ISCreateGeneral(PROTEUS_COMM_WORLD,elementBoundaries_subdomain_owned.size(),&elementBoundaryNumbering_new2old[0],PETSC_COPY_VALUES,&elementBoundaryNumberingIS_subdomain_new2old);
3711 IS elementBoundaryNumberingIS_global_new2old;
3712 ISAllGather(elementBoundaryNumberingIS_subdomain_new2old,&elementBoundaryNumberingIS_global_new2old);
3713 const PetscInt *elementBoundaryNumbering_global_new2old;
3715 ISGetIndices(elementBoundaryNumberingIS_global_new2old,&elementBoundaryNumbering_global_new2old);
3716 ierr =
enforceMemoryLimit(PROTEUS_COMM_WORLD, rank, max_rss_gb,
"Allocating elementBoudnary old2new/new2old");CHKERRABORT(PROTEUS_COMM_WORLD, ierr);
3719 elementBoundaryNumbering_global_old2new[elementBoundaryNumbering_global_new2old[ebN]] = ebN;
3721 ISRestoreIndices(elementBoundaryNumberingIS_global_new2old,&elementBoundaryNumbering_global_new2old);
3722 ISDestroy(&elementBoundaryNumberingIS_subdomain_new2old);
3723 ISDestroy(&elementBoundaryNumberingIS_global_new2old);
3724 PetscLogEventEnd(build_subdomains_faces_event,0,0,0,0);
3725 ierr =
enforceMemoryLimit(PROTEUS_COMM_WORLD, rank, max_rss_gb,
"Done allocating elementBoudnary old2new/new2old");CHKERRABORT(PROTEUS_COMM_WORLD, ierr);
3726 int build_subdomains_edges_event;
3727 PetscLogEventRegister(
"Subd edges",0,&build_subdomains_edges_event);
3728 PetscLogEventBegin(build_subdomains_edges_event,0,0,0,0);
3733 std::ifstream edgeFile(edgeFileName.c_str());
3735 if (!edgeFile.good())
3737 std::cerr<<
"cannot open Triangle edge file"
3738 <<edgeFileName<<std::endl;
3743 bool hasEdgeMarkers =
false;
3745 int ihasEdgeMarkers(0);
3747 assert(nEdges_global > 0);
3748 if (ihasEdgeMarkers > 0)
3750 hasEdgeMarkers =
true;
3754 set<int> edges_subdomain_owned;
3755 vector<set<int> > nodeEdgesStar(nNodes_subdomain_new[rank]);
3756 map<int,int> edgeMaterialTypesMap;
3757 map<int,vector<int> > elementEdgesMap;
3758 map<int,pair<int,int> > edgeNodesMap;
3759 set<int> supportedEdges;
3760 for (
int ied = 0; ied < nEdges_global; ied++)
3762 int ned,nn0,nn1;
int edId(0);
3764 if (ihasEdgeMarkers > 0)
3769 assert(0 <= ned && ned < nEdges_global && edgeFile.good());
3770 int nn0_new = nodeNumbering_global_old2new[nn0];
3771 if (nn0_new >= nodeOffsets_new[rank] && nn0_new < nodeOffsets_new[rank+1])
3773 nodeEdgesStar.at(nn0_new-nodeOffsets_new[rank]).insert(ned);
3774 supportedEdges.insert(ned);
3776 int nn1_new = nodeNumbering_global_old2new[nn1];
3777 if (nn1_new >= nodeOffsets_new[rank] && nn1_new < nodeOffsets_new[rank+1])
3779 nodeEdgesStar.at(nn1_new-nodeOffsets_new[rank]).insert(ned);
3780 supportedEdges.insert(ned);
3782 int nodes[2] = {nn0_new,nn1_new};
3785 if (edgeElementsMap.find(nodeTuple) != edgeElementsMap.end())
3787 if (nodeTuple.
nodes[0] >= nodeOffsets_new[rank] && nodeTuple.
nodes[0] < nodeOffsets_new[rank+1])
3788 edges_subdomain_owned.insert(ned);
3790 edgeNodesMap[ned].first = nodeTuple.
nodes[0];
3791 edgeNodesMap[ned].second = nodeTuple.
nodes[1];
3792 if (ihasEdgeMarkers > 0)
3793 edgeMaterialTypesMap[ned]=edId;
3794 for (set<pair<int,int> >::iterator elementp=edgeElementsMap[nodeTuple].begin();
3795 elementp != edgeElementsMap[nodeTuple].end();
3798 int eN = elementNumbering_global_old2new[elementp->first];
3799 if (elementEdgesMap.find(eN) != elementEdgesMap.end())
3801 elementEdgesMap[eN][elementp->second] = ned;
3805 std::vector<int>
init(3,-1);
3806 elementEdgesMap[eN] =
init;
3807 elementEdgesMap[eN][elementp->second] = ned;
3813 int nEdges_owned_subdomain=edges_subdomain_owned.size(),
3815 MPI_Allreduce(&nEdges_owned_subdomain,&nEdges_owned_new,1,MPI_INT,MPI_SUM,PROTEUS_COMM_WORLD);
3816 assert(nEdges_owned_new == nEdges_global);
3817 ierr =
enforceMemoryLimit(PROTEUS_COMM_WORLD, rank, max_rss_gb,
"Done reading edges");CHKERRABORT(PROTEUS_COMM_WORLD, ierr);
3823 for (map<
int,vector<int> >::iterator edgesp=elementEdgesMap.begin();
3824 edgesp!=elementEdgesMap.end();
3828 for (
int iv=0;iv<simplexDim;iv++)
3831 int nN_global = elementNodesArrayMap[elementNumbering_global_new2old[edgesp->first]][iv];
3832 if (nN_global >= nodeOffsets_new[rank] && nN_global < nodeOffsets_new[rank+1])
3835 for(
int ed=0;ed<3;ed++)
3837 nodeEdgesStar.at(nN_global-nodeOffsets_new[rank]).insert(edgesp->second[ed]);
3843 valarray<int> nodeEdgeOffsets(nNodes_subdomain_new[rank]+1);
3844 nodeEdgeOffsets[0] = 0;
3845 for (
int nN = 0; nN < nNodes_subdomain_new[rank]; nN++)
3846 nodeEdgeOffsets[nN+1] = nodeEdgeOffsets[nN]+nodeEdgesStar.at(nN).size();
3847 valarray<int> nodeEdgesArray(nodeEdgeOffsets[nNodes_subdomain_new[rank]]);
3848 for (
int nN=0,offset=0; nN < nNodes_subdomain_new[rank]; nN++)
3850 for (set<int>::iterator edN_star = nodeEdgesStar.at(nN).begin();
3851 edN_star != nodeEdgesStar.at(nN).end();
3852 edN_star++,offset++)
3854 nodeEdgesArray[offset] = *edN_star;
3859 valarray<int> nEdges_subdomain_new(size),
3860 edgeOffsets_new(size+1);
3861 for (
int sdN=0;sdN<size;sdN++)
3863 nEdges_subdomain_new[sdN] = edges_subdomain_owned.size();
3865 nEdges_subdomain_new[sdN] = 0;
3866 valarray<int> nEdges_subdomain_new_send=nEdges_subdomain_new;
3867 MPI_Allreduce(&nEdges_subdomain_new_send[0],&nEdges_subdomain_new[0],size,MPI_INT,MPI_SUM,PROTEUS_COMM_WORLD);
3871 edgeOffsets_new[0] = 0;
3872 for (
int sdN=0;sdN<size;sdN++)
3873 edgeOffsets_new[sdN+1] = edgeOffsets_new[sdN]+nEdges_subdomain_new[sdN];
3879 valarray<int> edgeNumbering_new2old(edges_subdomain_owned.size());
3880 set<int>::iterator edN_ownedp=edges_subdomain_owned.begin();
3881 for (
int edN=0;edN<int(edges_subdomain_owned.size());edN++,edN_ownedp++)
3883 edgeNumbering_new2old[edN] = *edN_ownedp;
3885 IS edgeNumberingIS_subdomain_new2old;
3886 ISCreateGeneral(PROTEUS_COMM_WORLD,edges_subdomain_owned.size(),&edgeNumbering_new2old[0],PETSC_COPY_VALUES,&edgeNumberingIS_subdomain_new2old);
3887 IS edgeNumberingIS_global_new2old;
3888 ISAllGather(edgeNumberingIS_subdomain_new2old,&edgeNumberingIS_global_new2old);
3889 const PetscInt *edgeNumbering_global_new2old;
3890 valarray<int> edgeNumbering_global_old2new(newMesh.
nEdges_global);
3891 ISGetIndices(edgeNumberingIS_global_new2old,&edgeNumbering_global_new2old);
3892 ierr =
enforceMemoryLimit(PROTEUS_COMM_WORLD, rank, max_rss_gb,
"Setting edgeNumering old2new/new2old");CHKERRABORT(PROTEUS_COMM_WORLD, ierr);
3895 edgeNumbering_global_old2new[edgeNumbering_global_new2old[edN]] = edN;
3897 ISRestoreIndices(edgeNumberingIS_global_new2old,&edgeNumbering_global_new2old);
3898 ISDestroy(&edgeNumberingIS_subdomain_new2old);
3899 ISDestroy(&edgeNumberingIS_global_new2old);
3900 ierr =
enforceMemoryLimit(PROTEUS_COMM_WORLD, rank, max_rss_gb,
"Done allocating edgeNumering old2new/new2old");CHKERRABORT(PROTEUS_COMM_WORLD, ierr);
3905 set<int> elements_overlap,nodes_overlap,elementBoundaries_overlap,edges_overlap;
3906 for (
int nN = 0; nN < nNodes_subdomain_new[rank]; nN++)
3909 for (
int offset = nodeStarOffsetsNew[nN];offset<nodeStarOffsetsNew[nN+1];offset++)
3911 int nN_point_global = nodeStarArrayNew[offset];
3912 bool offproc = nN_point_global < nodeOffsets_new[rank] || nN_point_global >= nodeOffsets_new[rank+1];
3914 nodes_overlap.insert(nN_point_global);
3917 for (
int eN_star_offset = nodeElementOffsets[nN];
3918 eN_star_offset < nodeElementOffsets[nN+1]; eN_star_offset++)
3920 int eN_star_old = nodeElementsArray[eN_star_offset];
3921 int eN_star_new = elementNumbering_global_old2new[eN_star_old];
3922 bool offproc = eN_star_new >= elementOffsets_new[rank+1] || eN_star_new < elementOffsets_new[rank];
3924 elements_overlap.insert(eN_star_new);
3927 for (
int ebN_star_offset = nodeElementBoundaryOffsets[nN];
3928 ebN_star_offset < nodeElementBoundaryOffsets[nN+1]; ebN_star_offset++)
3930 int ebN_star_old = nodeElementBoundariesArray[ebN_star_offset];
3931 int ebN_star_new = elementBoundaryNumbering_global_old2new[ebN_star_old];
3932 bool offproc = ebN_star_new >= elementBoundaryOffsets_new[rank+1] || ebN_star_new < elementBoundaryOffsets_new[rank];
3934 elementBoundaries_overlap.insert(ebN_star_new);
3937 for (
int edN_star_offset = nodeEdgeOffsets[nN];
3938 edN_star_offset < nodeEdgeOffsets[nN+1]; edN_star_offset++)
3940 int edN_star_old = nodeEdgesArray[edN_star_offset];
3941 int edN_star_new = edgeNumbering_global_old2new[edN_star_old];
3942 bool offproc = edN_star_new >= edgeOffsets_new[rank+1] || edN_star_new < edgeOffsets_new[rank];
3944 edges_overlap.insert(edN_star_new);
3947 elementNumbering_global_old2new.resize(0);
3948 MPI_Barrier(PROTEUS_COMM_WORLD);
3949 assert(edges_overlap.size() + nEdges_subdomain_new[rank] == edgeNodesMap.size());
3953 int nN_subdomain = nNodes_subdomain_new[rank];
3954 map<int,int> nodes_overlap_global2subdomainMap;
3955 for (set<int>::iterator nN_globalp=nodes_overlap.begin();nN_globalp != nodes_overlap.end(); nN_globalp++,nN_subdomain++)
3956 nodes_overlap_global2subdomainMap[*nN_globalp] = nN_subdomain;
3975 PetscLogEventEnd(build_subdomains_edges_event,0,0,0,0);
3976 int build_subdomains_renumber_event;
3977 PetscLogEventRegister(
"Subd's renumber",0,&build_subdomains_renumber_event);
3978 PetscLogEventBegin(build_subdomains_renumber_event,0,0,0,0);
3983 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;
3984 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;
4002 map<int,int> nodeNumbering_global2subdomainMap;
4003 map<int,int> elementBoundaryNumbering_global2subdomainMap;
4004 map<int,int> edgeNumbering_global2subdomainMap;
4010 for (
int iv = 0; iv < nNodes_global; iv++)
4012 int nv;
double x,y;
int nodeId(0);
4014 if (hasVertexMarkers > 0)
4015 vertexFile >> nodeId;
4017 assert(0 <= nv && nv < nNodes_global && vertexFile.good());
4018 int nN_global_new = nodeNumbering_global_old2new[nv];
4020 if (nN_global_new >= nodeOffsets_new[rank] && nN_global_new < nodeOffsets_new[rank+1])
4022 int nv_subdomain_new = nN_global_new - nodeOffsets_new[rank];
4023 nodeNumbering_subdomain2global[nv_subdomain_new] = nN_global_new;
4024 nodeNumbering_global2subdomainMap[nN_global_new] = nv_subdomain_new;
4028 if (hasVertexMarkers > 0)
4032 if (nodes_overlap.count(nN_global_new) == 1)
4034 int nv_subdomain_new = nodes_overlap_global2subdomainMap[nN_global_new];
4035 nodeNumbering_subdomain2global[nv_subdomain_new] = nN_global_new;
4036 nodeNumbering_global2subdomainMap[nN_global_new] = nv_subdomain_new;
4040 if (hasVertexMarkers > 0)
4046 ISRestoreIndices(nodeNumberingIS_global_old2new,&nodeNumbering_global_old2new);
4047 ISDestroy(&nodePartitioningIS_new);
4048 ISDestroy(&nodeNumberingIS_subdomain_old2new);
4049 ISDestroy(&nodeNumberingIS_global_old2new);
4051 ierr =
enforceMemoryLimit(PROTEUS_COMM_WORLD, rank, max_rss_gb,
"Done reading vertices");CHKERRABORT(PROTEUS_COMM_WORLD, ierr);
4060 for (
int eN = 0; eN < nElements_subdomain_new[rank]; eN++)
4062 int eN_global_new = elementOffsets_new[rank] + eN;
4063 int eN_global_old = elementNumbering_global_new2old[eN_global_new];
4064 elementNumbering_subdomain2global[eN] = eN_global_new;
4068 int nN_global_new = elementNodesArrayMap[eN_global_old][nN];
4069 int nN_subdomain = nodeNumbering_global2subdomainMap[nN_global_new];
4076 set<int>::iterator eN_p = elements_overlap.begin();
4077 for (
int eN = nElements_subdomain_new[rank]; eN < nElements_subdomain_new[rank] + int(elements_overlap.size()); eN++,eN_p++)
4079 int eN_global_new = *eN_p;
4080 int eN_global_old = elementNumbering_global_new2old[eN_global_new];
4081 elementNumbering_subdomain2global[eN] = eN_global_new;
4085 int nN_global_new = elementNodesArrayMap[eN_global_old][nN];
4086 int nN_subdomain = nodeNumbering_global2subdomainMap[nN_global_new];
4090 ISRestoreIndices(elementNumberingIS_global_new2old,&elementNumbering_global_new2old);
4091 ISDestroy(&elementNumberingIS_subdomain_new2old);
4092 ISDestroy(&elementNumberingIS_global_new2old);
4098 for (
int ebN=0; ebN < nElementBoundaries_subdomain_new[rank]; ebN++)
4100 int ebN_global = ebN + elementBoundaryOffsets_new[rank];
4101 elementBoundaryNumbering_subdomain2global[ebN]=ebN_global;
4102 elementBoundaryNumbering_global2subdomainMap[ebN_global] = ebN;
4107 set<int>::iterator ebN_p = elementBoundaries_overlap.begin();
4108 for(
int ebN=nElementBoundaries_subdomain_new[rank];ebN < nElementBoundaries_subdomain_new[rank] + int(elementBoundaries_overlap.size()); ebN++,ebN_p++)
4110 int ebN_global = *ebN_p;
4111 elementBoundaryNumbering_subdomain2global[ebN] = ebN_global;
4112 elementBoundaryNumbering_global2subdomainMap[ebN_global] = ebN;
4121 for (
int eN=0;eN<nElements_subdomain_new[rank];eN++)
4123 int eN_global = eN+elementOffsets_new[rank];
4127 elementBoundaryNumbering_global2subdomainMap[elementBoundaryNumbering_global_old2new[elementBoundariesMap[eN_global][ebN]]];
4133 eN_p = elements_overlap.begin();
4134 for (
int eN = nElements_subdomain_new[rank]; eN < nElements_subdomain_new[rank] + int(elements_overlap.size()); eN++,eN_p++)
4136 int eN_global_new = *eN_p;
4140 elementBoundaryNumbering_global2subdomainMap[elementBoundaryNumbering_global_old2new[elementBoundariesMap[eN_global_new][ebN]]];
4148 for (
int edN=0; edN < nEdges_subdomain_new[rank]; edN++)
4150 int edN_global = edN + edgeOffsets_new[rank];
4151 edgeNumbering_subdomain2global[edN]=edN_global;
4152 edgeNumbering_global2subdomainMap[edN_global] = edN;
4157 set<int>::iterator edN_p = edges_overlap.begin();
4158 for(
int edN=nEdges_subdomain_new[rank];edN < nEdges_subdomain_new[rank] + int(edges_overlap.size()); edN++,edN_p++)
4160 int edN_global = *edN_p;
4161 edgeNumbering_subdomain2global[edN] = edN_global;
4162 edgeNumbering_global2subdomainMap[edN_global] = edN;
4170 for (map<
int,pair<int,int> >::iterator edgep=edgeNodesMap.begin();
4171 edgep!=edgeNodesMap.end();
4174 int edN_global_old = edgep->first;
4175 int edN_global_new = edgeNumbering_global_old2new[edN_global_old];
4176 assert(edgeNumbering_global2subdomainMap.find(edN_global_new) != edgeNumbering_global2subdomainMap.end());
4177 int edN_subdomain = edgeNumbering_global2subdomainMap[edN_global_new];
4181 edgeNumbering_global_old2new.resize(0);
4202 using namespace std;
4207 int> elementBoundaryIds;
4214 register int nodes[3];
4219 if(elementBoundaryElements.find(ebt) != elementBoundaryElements.end())
4221 elementBoundaryElements[ebt].right=eN;
4222 elementBoundaryElements[ebt].right_ebN_element=ebN;
4223 assert(elementBoundaryIds[ebt] == ebN_global);
4227 elementBoundaryElements.insert(elementBoundaryElements.end(),make_pair(ebt,
ElementNeighbors(eN,ebN)));
4228 elementBoundaryIds.insert(elementBoundaryIds.end(),make_pair(ebt,ebN_global));
4238 set<int> interiorElementBoundaries,exteriorElementBoundaries;
4249 eb != elementBoundaryElements.end();
4252 int ebN = elementBoundaryIds[eb->first];
4262 if(eb->second.right != -1)
4264 interiorElementBoundaries.insert(ebN);
4268 exteriorElementBoundaries.insert(ebN);
4270 if (eb->second.right != -1)
4280 for (set<int>::iterator ebN=interiorElementBoundaries.begin();ebN != interiorElementBoundaries.end(); ebN++,ebNI++)
4282 for (set<int>::iterator ebN=exteriorElementBoundaries.begin();ebN != exteriorElementBoundaries.end(); ebN++,ebNE++)
4284 set<NodeTuple<2> > edges;
4309 for (set<int>::iterator nN_star=nodeStar[nN].begin();nN_star!=nodeStar[nN].end();nN_star++,offset++)
4329 for (set<int>::iterator eN_star = nodeElementsStar[nN].begin(); eN_star != nodeElementsStar[nN].end();
4368 if (hasElementBoundaryMarkers)
4371 for (map<int,int>::iterator ebmp = elementBoundaryMaterialTypesMap.begin(); ebmp != elementBoundaryMaterialTypesMap.end();ebmp++)
4373 int ebN_global_new = elementBoundaryNumbering_global_old2new[ebmp->first];
4374 assert(elementBoundaryNumbering_global2subdomainMap.find(ebN_global_new) != elementBoundaryNumbering_global2subdomainMap.end());
4375 int ebN_subdomain = elementBoundaryNumbering_global2subdomainMap[ebN_global_new];
4378 if (!hasVertexMarkers)
4389 elementBoundaryNumbering_global_old2new.resize(0);
4390 ierr =
enforceMemoryLimit(PROTEUS_COMM_WORLD, rank, max_rss_gb,
"Done with material types");CHKERRABORT(PROTEUS_COMM_WORLD, ierr);
4391 PetscLogEventEnd(build_subdomains_renumber_event,0,0,0,0);
4392 int build_subdomains_cleanup_event;
4393 PetscLogEventRegister(
"Cleanup",0,&build_subdomains_cleanup_event);
4394 PetscLogEventBegin(build_subdomains_cleanup_event,0,0,0,0);
4408 for (
int sdN = 0; sdN < size+1; sdN++)
4439 H5Sclose(filespace);
4444 PetscLogEventEnd(build_subdomains_cleanup_event,0,0,0,0);
4446 PetscLogView(PETSC_VIEWER_STDOUT_WORLD);
4447 ierr =
enforceMemoryLimit(PROTEUS_COMM_WORLD, rank, max_rss_gb,
"Done with partitioning!");CHKERRABORT(PROTEUS_COMM_WORLD, ierr);
4453 using namespace std;
4456 ierr = MPI_Comm_size(PROTEUS_COMM_WORLD,&size);
4457 ierr = MPI_Comm_rank(PROTEUS_COMM_WORLD,&rank);
4491 valarray<int> elementOffsets_old(size+1);
4492 elementOffsets_old[0] = 0;
4493 for(
int sdN=0;sdN<size;sdN++)
4494 elementOffsets_old[sdN+1] = elementOffsets_old[sdN] +
4500 int nElements_subdomain = (elementOffsets_old[rank+1] - elementOffsets_old[rank]);
4501 PetscInt *elementNeighborsOffsets_subdomain,*elementNeighbors_subdomain,*weights_subdomain;
4502 PetscMalloc(
sizeof(PetscInt)*(nElements_subdomain+1),&elementNeighborsOffsets_subdomain);
4506 elementNeighborsOffsets_subdomain[0] = 0;
4507 for (
int eN=0,offset=0; eN < nElements_subdomain; eN++)
4509 int eN_global = elementOffsets_old[rank] + eN;
4510 int offsetStart=offset;
4514 if (eN_neighbor_global >= 0 )
4515 elementNeighbors_subdomain[offset++] = eN_neighbor_global;
4517 elementNeighborsOffsets_subdomain[eN+1]=offset;
4518 sort(&elementNeighbors_subdomain[offsetStart],&elementNeighbors_subdomain[offset]);
4519 int weight = (elementNeighborsOffsets_subdomain[eN+1] - elementNeighborsOffsets_subdomain[eN]);
4520 for (
int k=elementNeighborsOffsets_subdomain[eN];k<elementNeighborsOffsets_subdomain[eN+1];k++)
4521 weights_subdomain[k] = weight;
4530 ierr = MatCreateMPIAdj(PROTEUS_COMM_WORLD,
4531 nElements_subdomain,
4533 elementNeighborsOffsets_subdomain,
4534 elementNeighbors_subdomain,
4536 &petscAdjacency);CHKERRABORT(PROTEUS_COMM_WORLD, ierr);
4537 PetscFree(weights_subdomain);
4538 MatPartitioning petscPartition;
4539 MatPartitioningCreate(PROTEUS_COMM_WORLD,&petscPartition);
4540 MatPartitioningSetAdjacency(petscPartition,petscAdjacency);
4541 MatPartitioningSetFromOptions(petscPartition);
4544 IS elementPartitioningIS_new;
4545 MatPartitioningApply(petscPartition,&elementPartitioningIS_new);
4546 MatPartitioningDestroy(&petscPartition);
4547 MatDestroy(&petscAdjacency);
4579 valarray<int> nElements_subdomain_new(size);
4580 ISPartitioningCount(elementPartitioningIS_new,size,&nElements_subdomain_new[0]);
4583 valarray<int> elementOffsets_new(size+1);
4584 elementOffsets_new[0] = 0;
4585 for (
int sdN=0;sdN<size;sdN++)
4586 elementOffsets_new[sdN+1] = elementOffsets_new[sdN] + nElements_subdomain_new[sdN];
4589 IS elementNumberingIS_subdomain_old2new;
4590 ISPartitioningToNumbering(elementPartitioningIS_new,&elementNumberingIS_subdomain_old2new);
4597 IS elementNumberingIS_global_old2new;
4598 ISAllGather(elementNumberingIS_subdomain_old2new,&elementNumberingIS_global_old2new);
4599 const PetscInt *elementNumbering_global_old2new;
4600 ISGetIndices(elementNumberingIS_global_old2new,&elementNumbering_global_old2new);
4603 elementNumbering_global_new2old[elementNumbering_global_old2new[eN]] = eN;
4641 elementNumbering_global_old2new[eN_ebN];
4653 MPI_Recv(nodeMask,PetscBTLength(mesh.
nNodes_global),MPI_CHAR,rank-1,0,PROTEUS_COMM_WORLD,&status);
4656 set<int> nodes_subdomain_owned;
4657 for(
int eN=elementOffsets_new[rank];eN<elementOffsets_new[rank+1];eN++)
4660 int nN_global = elementNodesArray_new[eN*mesh.
nNodes_element+nN];
4661 if (!PetscBTLookupSet(nodeMask,nN_global))
4662 nodes_subdomain_owned.insert(nN_global);
4666 MPI_Send(nodeMask,PetscBTLength(mesh.
nNodes_global),MPI_CHAR,rank+1,0,PROTEUS_COMM_WORLD);
4667 ierr = PetscBTDestroy(&nodeMask);
4669 cerr<<
"Error in PetscBTDestroy"<<endl;
4671 valarray<int> nNodes_subdomain_new(size),
4672 nodeOffsets_new(size+1);
4673 for (
int sdN=0;sdN<size;sdN++)
4675 nNodes_subdomain_new[sdN] = nodes_subdomain_owned.size();
4677 nNodes_subdomain_new[sdN] = 0;
4678 valarray<int> nNodes_subdomain_new_send=nNodes_subdomain_new;
4679 MPI_Allreduce(&nNodes_subdomain_new_send[0],&nNodes_subdomain_new[0],size,MPI_INT,MPI_SUM,PROTEUS_COMM_WORLD);
4680 nodeOffsets_new[0] = 0;
4681 for (
int sdN=0;sdN<size;sdN++)
4682 nodeOffsets_new[sdN+1] = nodeOffsets_new[sdN]+nNodes_subdomain_new[sdN];
4688 valarray<int> nodeNumbering_new2old(nodes_subdomain_owned.size());
4689 set<int>::iterator nN_ownedp=nodes_subdomain_owned.begin();
4690 for (
int nN=0;nN<int(nodes_subdomain_owned.size());nN++)
4692 nodeNumbering_new2old[nN] = *nN_ownedp++;
4694 IS nodeNumberingIS_new2old;
4695 ISCreateGeneral(PROTEUS_COMM_WORLD,nodes_subdomain_owned.size(),&nodeNumbering_new2old[0],PETSC_COPY_VALUES,&nodeNumberingIS_new2old);
4696 IS nodeNumberingIS_global_new2old;
4697 ISAllGather(nodeNumberingIS_new2old,&nodeNumberingIS_global_new2old);
4698 const PetscInt *nodeNumbering_global_new2old;
4699 valarray<int> nodeNumbering_old2new_global(mesh.
nNodes_global);
4700 ISGetIndices(nodeNumberingIS_global_new2old,&nodeNumbering_global_new2old);
4703 nodeNumbering_old2new_global[nodeNumbering_global_new2old[nN]] = nN;
4711 elementNodesArray_new[eN*mesh.
nNodes_element+nN] = nodeNumbering_old2new_global[nN_old];
4723 edgeNodesArray_new[i] = nodeNumbering_old2new_global[nN_old];
4729 nodeStarArray_new[i] = nodeNumbering_old2new_global[nN_old];
4735 int nN_new = nodeNumbering_old2new_global[nN];
4736 nodeArray_new[nN_new*3+0] = mesh.
nodeArray[nN*3+0];
4737 nodeArray_new[nN_new*3+1] = mesh.
nodeArray[nN*3+1];
4738 nodeArray_new[nN_new*3+2] = mesh.
nodeArray[nN*3+2];
4743 MPI_Status status_elementBoundaries;
4744 PetscBT elementBoundaryMask;
4748 MPI_Recv(elementBoundaryMask,PetscBTLength(mesh.
nElementBoundaries_global),MPI_CHAR,rank-1,0,PROTEUS_COMM_WORLD,&status_elementBoundaries);
4751 set<int> elementBoundaries_subdomain_owned;
4752 for(
int eN=elementOffsets_new[rank];eN<elementOffsets_new[rank+1];eN++)
4756 if(!PetscBTLookup(elementBoundaryMask,ebN_global))
4762 if (nodes_subdomain_owned.count(nN_global_old) > 0)
4770 PetscBTSet(elementBoundaryMask,ebN_global);
4771 elementBoundaries_subdomain_owned.insert(ebN_global);
4775 PetscBTSet(elementBoundaryMask,ebN_global);
4776 elementBoundaries_subdomain_owned.insert(ebN_global);
4784 ierr = PetscBTDestroy(&elementBoundaryMask);
4786 cerr<<
"Error in PetscBTDestroy for elementBoundaries"<<endl;
4788 valarray<int> nElementBoundaries_subdomain_new(size),
4789 elementBoundaryOffsets_new(size+1);
4790 for (
int sdN=0;sdN<size;sdN++)
4792 nElementBoundaries_subdomain_new[sdN] = elementBoundaries_subdomain_owned.size();
4794 nElementBoundaries_subdomain_new[sdN] = 0;
4795 valarray<int> nElementBoundaries_subdomain_new_send=nElementBoundaries_subdomain_new;
4796 MPI_Allreduce(&nElementBoundaries_subdomain_new_send[0],&nElementBoundaries_subdomain_new[0],size,MPI_INT,MPI_SUM,PROTEUS_COMM_WORLD);
4797 elementBoundaryOffsets_new[0] = 0;
4798 for (
int sdN=0;sdN<size;sdN++)
4799 elementBoundaryOffsets_new[sdN+1] = elementBoundaryOffsets_new[sdN]+nElementBoundaries_subdomain_new[sdN];
4804 valarray<int> elementBoundaryNumbering_new2old(elementBoundaries_subdomain_owned.size());
4805 set<int>::iterator ebN_ownedp=elementBoundaries_subdomain_owned.begin();
4806 for (
int ebN=0;ebN<int(elementBoundaries_subdomain_owned.size());ebN++)
4808 elementBoundaryNumbering_new2old[ebN] = *ebN_ownedp++;
4810 IS elementBoundaryNumberingIS_new2old;
4811 ISCreateGeneral(PROTEUS_COMM_WORLD,elementBoundaries_subdomain_owned.size(),&elementBoundaryNumbering_new2old[0],PETSC_COPY_VALUES,&elementBoundaryNumberingIS_new2old);
4812 IS elementBoundaryNumberingIS_global_new2old;
4813 ISAllGather(elementBoundaryNumberingIS_new2old,&elementBoundaryNumberingIS_global_new2old);
4814 const PetscInt *elementBoundaryNumbering_global_new2old;
4816 ISGetIndices(elementBoundaryNumberingIS_global_new2old,&elementBoundaryNumbering_global_new2old);
4819 elementBoundaryNumbering_old2new_global[elementBoundaryNumbering_global_new2old[ebN]] = ebN;
4833 int ebN_old=elementBoundaryNumbering_global_new2old[ebN];
4841 elementBoundaryElementsArray_new[ebN*2+0] = elementNumbering_global_old2new[eN_L_old];
4843 elementBoundaryElementsArray_new[ebN*2+1] = elementNumbering_global_old2new[eN_R_old];
4904 MPI_Status status_edges;
4909 MPI_Recv(edgesMask,PetscBTLength(mesh.
nEdges_global),MPI_CHAR,rank-1,0,PROTEUS_COMM_WORLD,&status_edges);
4912 map<NodeTuple<2>,
int> nodesEdgeMap_global;
4913 set<int> edges_subdomain_owned;
4917 nodes[0] = edgeNodesArray_new[2*ig];
4918 nodes[1] = edgeNodesArray_new[2*ig+1];
4920 nodesEdgeMap_global[et] = ig;
4922 for(
int eN=elementOffsets_new[rank];eN<elementOffsets_new[rank+1];eN++)
4931 if (nodesEdgeMap_global.find(et) != nodesEdgeMap_global.end())
4933 int edge_global = nodesEdgeMap_global[et];
4934 if(!PetscBTLookup(edgesMask,edge_global))
4936 PetscBTSet(edgesMask,edge_global);
4937 edges_subdomain_owned.insert(nodesEdgeMap_global[et]);
4943 MPI_Send(edgesMask,PetscBTLength(mesh.
nEdges_global),MPI_CHAR,rank+1,0,PROTEUS_COMM_WORLD);
4944 ierr = PetscBTDestroy(&edgesMask);
4946 cerr<<
"Error in PetscBTDestroy for edges"<<endl;
4948 valarray<int> nEdges_subdomain_new(size),
4949 edgeOffsets_new(size+1);
4951 for (
int sdN=0; sdN < size; sdN++)
4954 nEdges_subdomain_new[sdN] = edges_subdomain_owned.size();
4956 nEdges_subdomain_new[sdN] = 0;
4959 valarray<int> nEdges_subdomain_new_send=nEdges_subdomain_new;
4960 MPI_Allreduce(&nEdges_subdomain_new_send[0],&nEdges_subdomain_new[0],size,MPI_INT,MPI_SUM,PROTEUS_COMM_WORLD);
4961 edgeOffsets_new[0] = 0;
4962 for (
int sdN=0;sdN<size;sdN++)
4963 edgeOffsets_new[sdN+1] = edgeOffsets_new[sdN]+nEdges_subdomain_new[sdN];
4966 valarray<int> edgeNumbering_new2old(edges_subdomain_owned.size());
4967 set<int>::iterator edges_ownedp = edges_subdomain_owned.begin();
4968 for (
int i=0; i < int(edges_subdomain_owned.size());i++)
4969 edgeNumbering_new2old[i] = *edges_ownedp++;
4971 IS edgeNumberingIS_new2old;
4972 ISCreateGeneral(PROTEUS_COMM_WORLD,edges_subdomain_owned.size(),&edgeNumbering_new2old[0],PETSC_COPY_VALUES,&edgeNumberingIS_new2old);
4973 IS edgeNumberingIS_global_new2old;
4974 ISAllGather(edgeNumberingIS_new2old,&edgeNumberingIS_global_new2old);
4975 const PetscInt *edgeNumbering_global_new2old;
4978 valarray<int> edgeNumbering_old2new_global(mesh.
nEdges_global);
4979 ISGetIndices(edgeNumberingIS_global_new2old,&edgeNumbering_global_new2old);
4982 edgeNumbering_old2new_global[edgeNumbering_global_new2old[ig]] = ig;
4987 valarray<int> edgeNodesArray_newNodesAndEdges(2*mesh.
nEdges_global);
4988 map<NodeTuple<2>,
int > nodesEdgeMap_global_new;
4992 const int edge_old = edgeNumbering_global_new2old[ig];
4993 nodes[0] = edgeNodesArray_new[edge_old*2+0];
4994 nodes[1] = edgeNodesArray_new[edge_old*2+1];
4996 edgeNodesArray_newNodesAndEdges[ig*2+0] = edgeNodesArray_new[edge_old*2+0];
4997 edgeNodesArray_newNodesAndEdges[ig*2+1] = edgeNodesArray_new[edge_old*2+1];
4998 assert(nodesEdgeMap_global_new.find(et) == nodesEdgeMap_global_new.end());
4999 nodesEdgeMap_global_new[et] = ig;
5004 set<int> elements_overlap,nodes_overlap,elementBoundaries_overlap,edges_overlap;
5005 for(
int eN=elementOffsets_new[rank];eN<elementOffsets_new[rank+1];eN++)
5009 int nN_global = elementNodesArray_new[eN*mesh.
nNodes_element+nN];
5010 if (nN_global < nodeOffsets_new[rank] || nN_global >= nodeOffsets_new[rank+1])
5011 nodes_overlap.insert(nN_global);
5016 if (ebN_global < elementBoundaryOffsets_new[rank] || ebN_global >= elementBoundaryOffsets_new[rank+1])
5017 elementBoundaries_overlap.insert(ebN_global);
5026 const int nN0_global = elementNodesArray_new[eN*mesh.
nNodes_element+nN0];
5027 const int nN1_global = elementNodesArray_new[eN*mesh.
nNodes_element+nN1];
5028 if (nN0_global < nodeOffsets_new[rank] || nN0_global >= nodeOffsets_new[rank+1])
5029 assert(nodes_overlap.find(nN0_global) != nodes_overlap.end());
5030 if (nN1_global < nodeOffsets_new[rank] || nN1_global >= nodeOffsets_new[rank+1])
5031 assert(nodes_overlap.find(nN1_global) != nodes_overlap.end());
5032 if (nodesEdgeMap_global_new.find(et) != nodesEdgeMap_global_new.end())
5034 const int edge_global = nodesEdgeMap_global_new[et];
5035 if (edge_global < edgeOffsets_new[rank] || edge_global >= edgeOffsets_new[rank+1])
5036 edges_overlap.insert(edge_global);
5041 if (nElements_overlap > 0)
5064 for (set<int>::iterator nN=nodes_subdomain_owned.begin();nN != nodes_subdomain_owned.end();nN++)
5069 int eN_new = elementNumbering_global_old2new[eN];
5070 if (eN_new < elementOffsets_new[rank] or eN_new >= elementOffsets_new[rank+1])
5072 elements_overlap.insert(eN_new);
5073 for (
int nN_element=0;nN_element<mesh.
nNodes_element;nN_element++)
5075 int nN_global = elementNodesArray_new[eN_new*mesh.
nNodes_element+nN_element];
5076 if (nN_global < nodeOffsets_new[rank] or nN_global >= nodeOffsets_new[rank+1])
5077 nodes_overlap.insert(nN_global);
5082 if (ebN_global < elementBoundaryOffsets_new[rank] || ebN_global >= elementBoundaryOffsets_new[rank+1])
5083 elementBoundaries_overlap.insert(ebN_global);
5090 nodes[0] = elementNodesArray_new[eN_new*mesh.
nNodes_element+nN0];
5091 nodes[1] = elementNodesArray_new[eN_new*mesh.
nNodes_element+nN1];
5093 const int nN0_global = elementNodesArray_new[eN_new*mesh.
nNodes_element+nN0];
5094 const int nN1_global = elementNodesArray_new[eN_new*mesh.
nNodes_element+nN1];
5095 if (nodesEdgeMap_global_new.find(et) != nodesEdgeMap_global_new.end())
5097 const int edge_global = nodesEdgeMap_global_new[et];
5098 if (edge_global < edgeOffsets_new[rank] || edge_global >= edgeOffsets_new[rank+1])
5099 edges_overlap.insert(edge_global);
5106 for(
int eN=elementOffsets_new[rank];eN<elementOffsets_new[rank+1];eN++)
5111 (eN_ebN < elementOffsets_new[rank] || eN_ebN >= elementOffsets_new[rank+1]))
5113 elements_overlap.insert(eN_ebN);
5116 int nN_global = elementNodesArray_new[eN_ebN*mesh.
nNodes_element+nN];
5117 if (nN_global < nodeOffsets_new[rank] || nN_global >= nodeOffsets_new[rank+1])
5118 nodes_overlap.insert(nN_global);
5123 if (ebN_global < elementBoundaryOffsets_new[rank] || ebN_global >= elementBoundaryOffsets_new[rank+1])
5124 elementBoundaries_overlap.insert(ebN_global);
5131 nodes[0] = elementNodesArray_new[eN_ebN*mesh.
nNodes_element+nN0];
5132 nodes[1] = elementNodesArray_new[eN_ebN*mesh.
nNodes_element+nN1];
5134 const int nN0_global = elementNodesArray_new[eN_ebN*mesh.
nNodes_element+nN0];
5135 const int nN1_global = elementNodesArray_new[eN_ebN*mesh.
nNodes_element+nN1];
5136 bool foundEdge =
false;
5137 const int edge_global = nodesEdgeMap_global_new[et];
5138 if (edge_global < edgeOffsets_new[rank] || edge_global >= edgeOffsets_new[rank+1])
5139 edges_overlap.insert(edge_global);
5144 for (
int layer=1;layer<nElements_overlap;layer++)
5146 for (set<int>::iterator eN_p=elements_overlap.begin();eN_p != elements_overlap.end();eN_p++)
5148 int eN_global = *eN_p;
5153 (eN_ebN < elementOffsets_new[rank] || eN_ebN >= elementOffsets_new[rank+1]))
5155 elements_overlap.insert(eN_ebN);
5158 int nN_global = elementNodesArray_new[eN_ebN*mesh.
nNodes_element+nN];
5159 if (nN_global < nodeOffsets_new[rank] || nN_global >= nodeOffsets_new[rank+1])
5160 nodes_overlap.insert(nN_global);
5165 if (ebN_global < elementBoundaryOffsets_new[rank] || ebN_global >= elementBoundaryOffsets_new[rank+1])
5166 elementBoundaries_overlap.insert(ebN_global);
5173 nodes[0] = elementNodesArray_new[eN_ebN*mesh.
nNodes_element+nN0];
5174 nodes[1] = elementNodesArray_new[eN_ebN*mesh.
nNodes_element+nN1];
5176 const int nN0_global = elementNodesArray_new[eN_ebN*mesh.
nNodes_element+nN0];
5177 const int nN1_global = elementNodesArray_new[eN_ebN*mesh.
nNodes_element+nN1];
5178 bool foundEdge =
false;
5179 if (nodesEdgeMap_global_new.find(et) != nodesEdgeMap_global_new.end())
5181 const int edge_global = nodesEdgeMap_global_new[et];
5182 if (edge_global < edgeOffsets_new[rank] || edge_global >= edgeOffsets_new[rank+1])
5183 edges_overlap.insert(edge_global);
5207 map<int,int> nodeNumbering_global2subdomain;
5210 for(
int nN=0;nN<nNodes_subdomain_new[rank];nN++)
5212 int nN_global = nN + nodeOffsets_new[rank];
5213 nodeNumbering_subdomain2global[nN] = nN_global;
5214 nodeNumbering_global2subdomain[nN_global] = nN;
5222 set<int>::iterator nN_p=nodes_overlap.begin();
5223 for(
int nN=nNodes_subdomain_new[rank];nN < nNodes_subdomain_new[rank] + int(nodes_overlap.size()); nN++)
5225 int nN_global = *nN_p++;
5226 nodeNumbering_subdomain2global[nN] = nN_global;
5227 nodeNumbering_global2subdomain[nN_global] = nN;
5236 map<int,int> elementBoundaryNumbering_global2subdomain;
5237 for (
int ebN=0; ebN < nElementBoundaries_subdomain_new[rank]; ebN++)
5239 int ebN_global = ebN + elementBoundaryOffsets_new[rank];
5240 elementBoundaryNumbering_subdomain2global[ebN]=ebN_global;
5241 elementBoundaryNumbering_global2subdomain[ebN_global] = ebN;
5243 set<int>::iterator ebN_p = elementBoundaries_overlap.begin();
5244 for(
int ebN=nElementBoundaries_subdomain_new[rank];ebN < nElementBoundaries_subdomain_new[rank] + int(elementBoundaries_overlap.size()); ebN++)
5246 int ebN_global = *ebN_p++;
5247 elementBoundaryNumbering_subdomain2global[ebN] = ebN_global;
5248 elementBoundaryNumbering_global2subdomain[ebN_global] = ebN;
5258 for (
int eN=0;eN<nElements_subdomain_new[rank];eN++)
5260 int eN_global = eN+elementOffsets_new[rank];
5261 elementNumbering_subdomain2global[eN] = eN_global;
5266 nodeNumbering_global2subdomain[elementNodesArray_new[eN_global*mesh.
nNodes_element + nN]];
5272 set<int>::iterator eN_p=elements_overlap.begin();
5273 for(
int eN=nElements_subdomain_new[rank];eN < nElements_subdomain_new[rank]+int(elements_overlap.size());eN++)
5275 int eN_global = *eN_p++;
5277 elementNumbering_subdomain2global[eN] = eN_global;
5281 nodeNumbering_global2subdomain[elementNodesArray_new[eN_global*mesh.
nNodes_element + nN]];
5292 for (
int i=0; i < nEdges_subdomain_new[rank]; i++)
5294 const int ig = i+edgeOffsets_new[rank];
5295 const int nN0_global = edgeNodesArray_newNodesAndEdges[ig*2+0];
5296 const int nN1_global = edgeNodesArray_newNodesAndEdges[ig*2+1];
5298 assert(nodeNumbering_global2subdomain.find(nN0_global) != nodeNumbering_global2subdomain.end());
5299 assert(nodeNumbering_global2subdomain.find(nN1_global) != nodeNumbering_global2subdomain.end());
5300 const int nN0_subdomain = nodeNumbering_global2subdomain[nN0_global];
5301 const int nN1_subdomain = nodeNumbering_global2subdomain[nN1_global];
5304 edgeNumbering_subdomain2global[i] = ig;
5306 set<int>::iterator edge_p = edges_overlap.begin();
5307 for (
int i=nEdges_subdomain_new[rank]; i < nEdges_subdomain_new[rank] + int(edges_overlap.size()); i++)
5309 const int ig =*edge_p++;
5310 const int nN0_global = edgeNodesArray_newNodesAndEdges[ig*2+0];
5311 const int nN1_global = edgeNodesArray_newNodesAndEdges[ig*2+1];
5313 const int nN0_subdomain = nodeNumbering_global2subdomain[nN0_global];
5314 const int nN1_subdomain = nodeNumbering_global2subdomain[nN1_global];
5317 edgeNumbering_subdomain2global[i] = ig;
5378 int eN_global_new = elementNumbering_subdomain2global[eN];
5379 int eN_global_old = elementNumbering_global_new2old[eN_global_new];
5395 for (
int sdN=0;sdN<size+1;sdN++)
5462 ISRestoreIndices(elementNumberingIS_global_old2new,&elementNumbering_global_old2new);
5464 ISDestroy(&elementPartitioningIS_new);
5465 ISDestroy(&elementNumberingIS_subdomain_old2new);
5466 ISDestroy(&elementNumberingIS_global_old2new);
5468 ISRestoreIndices(nodeNumberingIS_global_new2old,&nodeNumbering_global_new2old);
5470 ISDestroy(&nodeNumberingIS_new2old);
5471 ISDestroy(&nodeNumberingIS_global_new2old);
5473 ISRestoreIndices(elementBoundaryNumberingIS_global_new2old,&elementBoundaryNumbering_global_new2old);
5475 ISDestroy(&elementBoundaryNumberingIS_new2old);
5476 ISDestroy(&elementBoundaryNumberingIS_global_new2old);
5478 ISRestoreIndices(edgeNumberingIS_global_new2old,&edgeNumbering_global_new2old);
5480 ISDestroy(&edgeNumberingIS_new2old);
5481 ISDestroy(&edgeNumberingIS_global_new2old);
5487 const int *elementOffsets_subdomain_owned,
5488 const int *nodeOffsets_subdomain_owned,
5489 const int *elementNumbering_subdomain2global,
5490 const int *nodeNumbering_subdomain2global,
5491 int& nDOF_all_processes,
5492 int& nDOF_subdomain,
5493 int& max_dof_neighbors,
5494 int *offsets_subdomain_owned,
5495 int * subdomain_l2g,
5496 int* subdomain2global,
5497 double * lagrangeNodesArray)
5499 using namespace std;
5502 ierr = MPI_Comm_size(PROTEUS_COMM_WORLD,&size);
5503 ierr = MPI_Comm_rank(PROTEUS_COMM_WORLD,&rank);
5510 assert(elementOffsets_subdomain_owned);
5511 assert(nodeOffsets_subdomain_owned);
5512 assert(elementNumbering_subdomain2global);
5513 assert(nodeNumbering_subdomain2global);
5516 const int nNodes_owned = nodeOffsets_subdomain_owned[rank+1]-nodeOffsets_subdomain_owned[rank];
5517 const int nElements_owned = elementOffsets_subdomain_owned[rank+1]-elementOffsets_subdomain_owned[rank];
5523 valarray<int> quadraticNumbering_new2old(nNodes_owned + nElements_owned);
5524 for (
int nN = 0; nN < nNodes_owned; nN++)
5526 int dof_global = nodeNumbering_subdomain2global[nN];
5527 quadraticNumbering_new2old[nN] = dof_global;
5529 for (
int eN = 0; eN < nElements_owned; eN++)
5531 int dof_global = mesh.
nNodes_global + elementNumbering_subdomain2global[eN];
5532 quadraticNumbering_new2old[nNodes_owned + eN] = dof_global;
5536 IS quadraticNumberingIS_new2old;
5537 ISCreateGeneral(PROTEUS_COMM_WORLD,nNodes_owned + nElements_owned,&quadraticNumbering_new2old[0],PETSC_COPY_VALUES,&quadraticNumberingIS_new2old);
5538 IS quadraticNumberingIS_global_new2old;
5539 ISAllGather(quadraticNumberingIS_new2old,&quadraticNumberingIS_global_new2old);
5541 const PetscInt *quadraticNumbering_global_new2old;
5543 ISGetIndices(quadraticNumberingIS_global_new2old,&quadraticNumbering_global_new2old);
5545 quadraticNumbering_old2new_global[quadraticNumbering_global_new2old[
id]] =
id;
5553 assert(offsets_subdomain_owned);
5554 assert(subdomain2global);
5555 assert(subdomain_l2g);
5556 assert(lagrangeNodesArray);
5561 for (
int sdN=0; sdN < size+1; sdN++)
5562 offsets_subdomain_owned[sdN] = nodeOffsets_subdomain_owned[sdN]+elementOffsets_subdomain_owned[sdN];
5566 int localOffset = 0;
5567 for (
int nN = 0; nN < nNodes_owned; nN++)
5569 int dof_global_old = nodeNumbering_subdomain2global[nN];
5570 int dof_global_new = quadraticNumbering_old2new_global[dof_global_old];
5571 subdomain2global[localOffset + nN] = dof_global_new;
5575 localOffset += nNodes_owned;
5576 for (
int eN = 0; eN < nElements_owned; eN++)
5578 int dof_global_old = mesh.
nNodes_global + elementNumbering_subdomain2global[eN];
5579 int dof_global_new = quadraticNumbering_old2new_global[dof_global_old];
5580 subdomain2global[localOffset + eN] = dof_global_new;
5584 localOffset += nElements_owned;
5587 int dof_global_old = nodeNumbering_subdomain2global[nN];
5588 int dof_global_new = quadraticNumbering_old2new_global[dof_global_old];
5589 subdomain2global[localOffset + nN -nNodes_owned] = dof_global_new;
5596 int dof_global_old = mesh.
nNodes_global + elementNumbering_subdomain2global[eN];
5597 int dof_global_new = quadraticNumbering_old2new_global[dof_global_old];
5598 subdomain2global[localOffset + eN-nElements_owned] = dof_global_new;
5603 const int nDOF_element = 3;
5604 const int ghostNodeOffset = nNodes_owned + nElements_owned;
5608 int nN_global_subdomain[2];
5614 if (nN_global_subdomain[nN] < nNodes_owned)
5615 subdomain_l2g[eN*nDOF_element+nN] = nN_global_subdomain[nN];
5617 subdomain_l2g[eN*nDOF_element+nN] = nN_global_subdomain[nN]-nNodes_owned + ghostNodeOffset;
5619 for (
int eI=0; eI < 3; eI++)
5620 lagrangeNodesArray[subdomain_l2g[eN*nDOF_element+nN]*3+eI] = mesh.
subdomainp->
nodeArray[3*nN_global_subdomain[nN]+eI];
5622 if (eN < nElements_owned)
5623 subdomain_l2g[eN*nDOF_element+2] = nNodes_owned + eN;
5625 subdomain_l2g[eN*nDOF_element+2] = eN - nElements_owned + ghostElementOffset;
5627 for (
int eI=0; eI < 3; eI++)
5631 ISRestoreIndices(quadraticNumberingIS_global_new2old,&quadraticNumbering_global_new2old);
5632 ISDestroy(&quadraticNumberingIS_new2old);
5633 ISDestroy(&quadraticNumberingIS_global_new2old);
5639 const int *elementBoundaryOffsets_subdomain_owned,
5640 const int *nodeOffsets_subdomain_owned,
5641 const int *elementBoundaryNumbering_subdomain2global,
5642 const int *nodeNumbering_subdomain2global,
5643 int& nDOF_all_processes,
5644 int& nDOF_subdomain,
5645 int& max_dof_neighbors,
5646 int *offsets_subdomain_owned,
5648 int *subdomain2global,
5649 double * lagrangeNodesArray)
5651 using namespace std;
5654 ierr = MPI_Comm_size(PROTEUS_COMM_WORLD,&size);
5655 ierr = MPI_Comm_rank(PROTEUS_COMM_WORLD,&rank);
5661 assert(elementBoundaryOffsets_subdomain_owned);
5662 assert(nodeOffsets_subdomain_owned);
5663 assert(elementBoundaryNumbering_subdomain2global);
5664 assert(nodeNumbering_subdomain2global);
5667 const int nNodes_owned = nodeOffsets_subdomain_owned[rank+1]-nodeOffsets_subdomain_owned[rank];
5668 const int nElementBoundaries_owned = elementBoundaryOffsets_subdomain_owned[rank+1]-elementBoundaryOffsets_subdomain_owned[rank];
5669 const int nDOFs_owned = nNodes_owned + nElementBoundaries_owned;
5675 valarray<int> quadraticNumbering_new2old(nDOFs_owned);
5676 for (
int nN = 0; nN < nNodes_owned; nN++)
5678 int dof_global = nodeNumbering_subdomain2global[nN];
5679 quadraticNumbering_new2old[nN] = dof_global;
5681 for (
int ebN = 0; ebN < nElementBoundaries_owned; ebN++)
5683 int dof_global = mesh.
nNodes_global + elementBoundaryNumbering_subdomain2global[ebN];
5684 quadraticNumbering_new2old[nNodes_owned + ebN] = dof_global;
5688 IS quadraticNumberingIS_new2old;
5689 ISCreateGeneral(PROTEUS_COMM_WORLD,nDOFs_owned,&quadraticNumbering_new2old[0],PETSC_COPY_VALUES,&quadraticNumberingIS_new2old);
5690 IS quadraticNumberingIS_global_new2old;
5691 ISAllGather(quadraticNumberingIS_new2old,&quadraticNumberingIS_global_new2old);
5693 const PetscInt *quadraticNumbering_global_new2old;
5694 valarray<int> quadraticNumbering_old2new_global(nDOFs_global);
5695 ISGetIndices(quadraticNumberingIS_global_new2old,&quadraticNumbering_global_new2old);
5696 for (
int id = 0;
id < nDOFs_global;
id++)
5697 quadraticNumbering_old2new_global[quadraticNumbering_global_new2old[
id]] =
id;
5705 assert(offsets_subdomain_owned);
5706 assert(subdomain2global);
5707 assert(subdomain_l2g);
5708 assert(lagrangeNodesArray);
5709 for (
int sdN=0; sdN < size+1; sdN++)
5710 offsets_subdomain_owned[sdN] = nodeOffsets_subdomain_owned[sdN]+elementBoundaryOffsets_subdomain_owned[sdN];
5717 int localOffset = 0;
5718 for (
int nN = 0; nN < nNodes_owned; nN++)
5720 int dof_global_old = nodeNumbering_subdomain2global[nN];
5721 int dof_global_new = quadraticNumbering_old2new_global[dof_global_old];
5722 subdomain2global[localOffset + nN] = dof_global_new;
5726 localOffset += nNodes_owned;
5727 for (
int ebN = 0; ebN < nElementBoundaries_owned; ebN++)
5729 int dof_global_old = mesh.
nNodes_global + elementBoundaryNumbering_subdomain2global[ebN];
5730 int dof_global_new = quadraticNumbering_old2new_global[dof_global_old];
5731 subdomain2global[localOffset + ebN] = dof_global_new;
5736 localOffset += nElementBoundaries_owned;
5739 int dof_global_old = nodeNumbering_subdomain2global[nN];
5740 int dof_global_new = quadraticNumbering_old2new_global[dof_global_old];
5741 subdomain2global[localOffset + nN -nNodes_owned] = dof_global_new;
5748 int dof_global_old = mesh.
nNodes_global + elementBoundaryNumbering_subdomain2global[ebN];
5749 int dof_global_new = quadraticNumbering_old2new_global[dof_global_old];
5750 subdomain2global[localOffset + ebN-nElementBoundaries_owned] = dof_global_new;
5755 const int nDOF_element = 6;
5756 const int ghostNodeOffset = nNodes_owned + nElementBoundaries_owned;
5760 int nN_global_subdomain[3];
5766 if (nN_global_subdomain[nN] < nNodes_owned)
5767 subdomain_l2g[eN*nDOF_element+nN] = nN_global_subdomain[nN];
5769 subdomain_l2g[eN*nDOF_element+nN] = nN_global_subdomain[nN]-nNodes_owned + ghostNodeOffset;
5771 for (
int eI=0; eI < 3; eI++)
5772 lagrangeNodesArray[subdomain_l2g[eN*nDOF_element+nN]*3+eI] = mesh.
subdomainp->
nodeArray[3*nN_global_subdomain[nN]+eI];
5779 if (ebN_global < nElementBoundaries_owned)
5780 subdomain_l2g[eN*nDOF_element+3+ebN] = nNodes_owned + ebN_global;
5782 subdomain_l2g[eN*nDOF_element+3+ebN] = ebN_global - nElementBoundaries_owned + ghostElementBoundaryOffset;
5784 for (
int eI=0; eI < 3; eI++)
5785 lagrangeNodesArray[subdomain_l2g[eN*nDOF_element+3+ebN]*3+eI] = 0.5*(mesh.
subdomainp->
nodeArray[3*nN_global_subdomain[(ebN+0)%3]+eI]+mesh.
subdomainp->
nodeArray[3*nN_global_subdomain[(ebN+1)%3]+eI]);
5789 ISRestoreIndices(quadraticNumberingIS_global_new2old,&quadraticNumbering_global_new2old);
5790 ISDestroy(&quadraticNumberingIS_new2old);
5791 ISDestroy(&quadraticNumberingIS_global_new2old);
5798 const int *edgeOffsets_subdomain_owned,
5799 const int *nodeOffsets_subdomain_owned,
5800 const int *edgeNumbering_subdomain2global,
5801 const int *nodeNumbering_subdomain2global,
5802 int& nDOF_all_processes,
5803 int& nDOF_subdomain,
5804 int& max_dof_neighbors,
5805 int *offsets_subdomain_owned,
5807 int *subdomain2global,
5808 double * lagrangeNodesArray)
5810 using namespace std;
5813 ierr = MPI_Comm_size(PROTEUS_COMM_WORLD,&size);
5814 ierr = MPI_Comm_rank(PROTEUS_COMM_WORLD,&rank);
5820 assert(edgeOffsets_subdomain_owned);
5821 assert(nodeOffsets_subdomain_owned);
5822 assert(edgeNumbering_subdomain2global);
5823 assert(nodeNumbering_subdomain2global);
5826 const int nNodes_owned = nodeOffsets_subdomain_owned[rank+1]-nodeOffsets_subdomain_owned[rank];
5827 const int nEdges_owned = edgeOffsets_subdomain_owned[rank+1]-edgeOffsets_subdomain_owned[rank];
5828 const int nDOFs_owned = nNodes_owned + nEdges_owned;
5834 valarray<int> quadraticNumbering_new2old(nDOFs_owned);
5835 for (
int nN = 0; nN < nNodes_owned; nN++)
5837 int dof_global = nodeNumbering_subdomain2global[nN];
5838 quadraticNumbering_new2old[nN] = dof_global;
5840 for (
int i = 0; i < nEdges_owned; i++)
5842 int dof_global = mesh.
nNodes_global + edgeNumbering_subdomain2global[i];
5843 quadraticNumbering_new2old[nNodes_owned + i] = dof_global;
5847 IS quadraticNumberingIS_new2old;
5848 ISCreateGeneral(PROTEUS_COMM_WORLD,nDOFs_owned,&quadraticNumbering_new2old[0],PETSC_COPY_VALUES,&quadraticNumberingIS_new2old);
5849 IS quadraticNumberingIS_global_new2old;
5850 ISAllGather(quadraticNumberingIS_new2old,&quadraticNumberingIS_global_new2old);
5852 const PetscInt *quadraticNumbering_global_new2old;
5853 valarray<int> quadraticNumbering_old2new_global(nDOFs_global);
5854 ISGetIndices(quadraticNumberingIS_global_new2old,&quadraticNumbering_global_new2old);
5855 for (
int id = 0;
id < nDOFs_global;
id++)
5856 quadraticNumbering_old2new_global[quadraticNumbering_global_new2old[
id]] =
id;
5867 assert(offsets_subdomain_owned);
5868 assert(subdomain2global);
5869 assert(subdomain_l2g);
5870 assert(lagrangeNodesArray);
5871 for (
int sdN=0; sdN < size+1; sdN++)
5872 offsets_subdomain_owned[sdN] = nodeOffsets_subdomain_owned[sdN]+edgeOffsets_subdomain_owned[sdN];
5879 int localOffset = 0;
5880 for (
int nN = 0; nN < nNodes_owned; nN++)
5882 int dof_global_old = nodeNumbering_subdomain2global[nN];
5883 int dof_global_new = quadraticNumbering_old2new_global[dof_global_old];
5884 subdomain2global[localOffset + nN] = dof_global_new;
5889 localOffset += nNodes_owned;
5890 for (
int i = 0; i < nEdges_owned; i++)
5892 int dof_global_old = mesh.
nNodes_global + edgeNumbering_subdomain2global[i];
5893 int dof_global_new = quadraticNumbering_old2new_global[dof_global_old];
5894 subdomain2global[localOffset + i] = dof_global_new;
5900 localOffset += nEdges_owned;
5903 int dof_global_old = nodeNumbering_subdomain2global[nN];
5904 int dof_global_new = quadraticNumbering_old2new_global[dof_global_old];
5905 subdomain2global[localOffset + nN -nNodes_owned] = dof_global_new;
5913 int dof_global_old = mesh.
nNodes_global + edgeNumbering_subdomain2global[i];
5914 int dof_global_new = quadraticNumbering_old2new_global[dof_global_old];
5915 subdomain2global[localOffset + i-nEdges_owned] = dof_global_new;
5921 const int nDOF_element = 10;
5922 const int ghostNodeOffset = nNodes_owned + nEdges_owned;
5925 map<NodeTuple<2>,
int> nodesEdgeMap_subdomain;
5934 nodesEdgeMap_subdomain[et] = i;
5938 int local_offset = 0;
5944 if (nN_global_subdomain < nNodes_owned)
5946 subdomain_l2g[eN*nDOF_element + local_offset + nN] = nN_global_subdomain;
5950 subdomain_l2g[eN*nDOF_element + local_offset + nN] = ghostNodeOffset + nN_global_subdomain - nNodes_owned;
5960 for (
int eI=0; eI < 3; eI++)
5961 lagrangeNodesArray[subdomain_l2g[eN*nDOF_element + local_offset +nN]*3+eI] = mesh.
subdomainp->
nodeArray[nN_global_subdomain*3+eI];
5971 int edge_subdomain = -1,edge_global=-1;
5974 nodes[0] = nN_global_subdomain;
5975 nodes[1] = nN_neig_global_subdomain;
5977 edge_subdomain = nodesEdgeMap_subdomain[et];
5978 edge_global = edgeNumbering_subdomain2global[edge_subdomain];
5979 assert(edge_subdomain >= 0 && edge_global >= 0);
5982 if (edge_subdomain < nEdges_owned)
5983 subdomain_l2g[eN*nDOF_element + local_offset + nN] = nNodes_owned + edge_subdomain;
5985 subdomain_l2g[eN*nDOF_element + local_offset + nN] = ghostEdgeOffset + edge_subdomain - nEdges_owned;
5994 for (
int eI=0; eI < 3; eI++)
5995 lagrangeNodesArray[subdomain_l2g[eN*nDOF_element+local_offset+nN]*3+eI] = 0.5*(mesh.
subdomainp->
nodeArray[nN_global_subdomain*3+eI]+
6006 int edge_subdomain = -1, edge_global = -1;
6009 nodes[0] = nN_global_subdomain;
6010 nodes[1] = nN_neig_global_subdomain;
6012 edge_subdomain = nodesEdgeMap_subdomain[et];
6013 edge_global = edgeNumbering_subdomain2global[edge_subdomain];
6014 assert(edge_subdomain >= 0 && edge_global >= 0);
6016 if (edge_subdomain < nEdges_owned)
6017 subdomain_l2g[eN*nDOF_element + local_offset + nN] = nNodes_owned + edge_subdomain;
6019 subdomain_l2g[eN*nDOF_element + local_offset + nN] = ghostEdgeOffset + edge_subdomain - nEdges_owned;
6028 for (
int eI=0; eI < 3; eI++)
6029 lagrangeNodesArray[subdomain_l2g[eN*nDOF_element+local_offset+nN]*3+eI] = 0.5*(mesh.
subdomainp->
nodeArray[nN_global_subdomain*3+eI]+
6040 int edge_subdomain = -1, edge_global = -1;
6043 nodes[0] = nN_global_subdomain;
6044 nodes[1] = nN_neig_global_subdomain;
6046 edge_subdomain = nodesEdgeMap_subdomain[et];
6047 edge_global = edgeNumbering_subdomain2global[edge_subdomain];
6048 assert(edge_subdomain >= 0 && edge_global >= 0);
6050 if (edge_subdomain < nEdges_owned)
6051 subdomain_l2g[eN*nDOF_element + local_offset + nN] = nNodes_owned + edge_subdomain;
6053 subdomain_l2g[eN*nDOF_element + local_offset + nN] = ghostEdgeOffset + edge_subdomain - nEdges_owned;
6062 for (
int eI=0; eI < 3; eI++)
6063 lagrangeNodesArray[subdomain_l2g[eN*nDOF_element+local_offset+nN]*3+eI] = 0.5*(mesh.
subdomainp->
nodeArray[nN_global_subdomain*3+eI]+
6069 ISRestoreIndices(quadraticNumberingIS_global_new2old,&quadraticNumbering_global_new2old);
6070 ISDestroy(&quadraticNumberingIS_new2old);
6071 ISDestroy(&quadraticNumberingIS_global_new2old);
6077 const int *edgeOffsets_subdomain_owned,
6078 const int *nodeOffsets_subdomain_owned,
6079 const int *edgeNumbering_subdomain2global,
6080 const int *nodeNumbering_subdomain2global,
6081 int& nDOF_all_processes,
6082 int& nDOF_subdomain,
6083 int& max_dof_neighbors,
6084 int *offsets_subdomain_owned,
6086 int *subdomain2global,
6087 double * lagrangeNodesArray)
6089 using namespace std;
6092 ierr = MPI_Comm_size(PROTEUS_COMM_WORLD,&size);
6093 ierr = MPI_Comm_rank(PROTEUS_COMM_WORLD,&rank);
6099 assert(edgeOffsets_subdomain_owned);
6100 assert(nodeOffsets_subdomain_owned);
6101 assert(edgeNumbering_subdomain2global);
6102 assert(nodeNumbering_subdomain2global);
6110 const int nNodes_owned = nodeOffsets_subdomain_owned[rank+1]-nodeOffsets_subdomain_owned[rank];
6111 const int nEdges_owned = edgeOffsets_subdomain_owned[rank+1]-edgeOffsets_subdomain_owned[rank];
6112 const int nBoundaries_owned = elementBoundaryOffsets_subdomain_owned[rank+1]-elementBoundaryOffsets_subdomain_owned[rank];
6113 const int nElements_owned = elementOffsets_subdomain_owned[rank+1]-elementOffsets_subdomain_owned[rank];
6115 const int nDOFs_owned = nNodes_owned + nEdges_owned + nBoundaries_owned + nElements_owned;
6121 valarray<int> quadraticNumbering_new2old(nDOFs_owned);
6122 for (
int nN = 0; nN < nNodes_owned; nN++)
6124 int dof_global = nodeNumbering_subdomain2global[nN];
6125 quadraticNumbering_new2old[nN] = dof_global;
6127 for (
int i = 0; i < nEdges_owned; i++)
6129 int dof_global = mesh.
nNodes_global + edgeNumbering_subdomain2global[i];
6130 quadraticNumbering_new2old[nNodes_owned + i] = dof_global;
6134 for (
int i = 0; i < nBoundaries_owned; i++)
6137 quadraticNumbering_new2old[nNodes_owned + nEdges_owned + i] = dof_global;
6139 for (
int i = 0; i < nElements_owned; i++)
6142 quadraticNumbering_new2old[nNodes_owned + nEdges_owned + nBoundaries_owned + i] = dof_global;
6148 IS quadraticNumberingIS_new2old;
6149 ISCreateGeneral(PROTEUS_COMM_WORLD,nDOFs_owned,&quadraticNumbering_new2old[0],PETSC_COPY_VALUES,&quadraticNumberingIS_new2old);
6150 IS quadraticNumberingIS_global_new2old;
6151 ISAllGather(quadraticNumberingIS_new2old,&quadraticNumberingIS_global_new2old);
6153 const PetscInt *quadraticNumbering_global_new2old;
6154 valarray<int> quadraticNumbering_old2new_global(nDOFs_global);
6155 ISGetIndices(quadraticNumberingIS_global_new2old,&quadraticNumbering_global_new2old);
6156 for (
int id = 0;
id < nDOFs_global;
id++)
6157 quadraticNumbering_old2new_global[quadraticNumbering_global_new2old[
id]] =
id;
6168 assert(offsets_subdomain_owned);
6169 assert(subdomain2global);
6170 assert(subdomain_l2g);
6171 assert(lagrangeNodesArray);
6172 for (
int sdN=0; sdN < size+1; sdN++)
6173 offsets_subdomain_owned[sdN] = nodeOffsets_subdomain_owned[sdN]+edgeOffsets_subdomain_owned[sdN]+elementBoundaryOffsets_subdomain_owned[sdN]+elementOffsets_subdomain_owned[sdN];
6180 int localOffset = 0;
6181 for (
int nN = 0; nN < nNodes_owned; nN++)
6183 int dof_global_old = nodeNumbering_subdomain2global[nN];
6184 int dof_global_new = quadraticNumbering_old2new_global[dof_global_old];
6185 subdomain2global[localOffset + nN] = dof_global_new;
6191 localOffset += nNodes_owned;
6192 for (
int i = 0; i < nEdges_owned; i++)
6194 int dof_global_old = mesh.
nNodes_global + edgeNumbering_subdomain2global[i];
6195 int dof_global_new = quadraticNumbering_old2new_global[dof_global_old];
6196 subdomain2global[localOffset + i] = dof_global_new;
6204 localOffset += nEdges_owned;
6205 for (
int i = 0; i < nBoundaries_owned; i++)
6208 int dof_global_new = quadraticNumbering_old2new_global[dof_global_old];
6209 subdomain2global[localOffset + i] = dof_global_new;
6216 localOffset += nBoundaries_owned;
6217 for (
int i = 0; i < nElements_owned; i++)
6220 int dof_global_new = quadraticNumbering_old2new_global[dof_global_old];
6221 subdomain2global[localOffset + i] = dof_global_new;
6228 localOffset += nElements_owned;
6231 int dof_global_old = nodeNumbering_subdomain2global[nN];
6232 int dof_global_new = quadraticNumbering_old2new_global[dof_global_old];
6233 subdomain2global[localOffset + nN - nNodes_owned] = dof_global_new;
6242 int dof_global_old = mesh.
nNodes_global + edgeNumbering_subdomain2global[i];
6243 int dof_global_new = quadraticNumbering_old2new_global[dof_global_old];
6244 subdomain2global[localOffset + i - nEdges_owned] = dof_global_new;
6254 int dof_global_new = quadraticNumbering_old2new_global[dof_global_old];
6255 subdomain2global[localOffset + i - nBoundaries_owned] = dof_global_new;
6265 int dof_global_new = quadraticNumbering_old2new_global[dof_global_old];
6266 subdomain2global[localOffset + i - nElements_owned] = dof_global_new;
6273 const int nDOF_element = 27;
6274 const int ghostNodeOffset = nNodes_owned + nEdges_owned + nBoundaries_owned + nElements_owned;
6279 int ledge[12][2] = {{0,1},{1,2},{2,3},{3,0},
6280 {0,4},{1,5},{2,6},{3,7},
6281 {4,5},{5,6},{6,7},{7,4}};
6283 int nEdges_element = 12;
6285 int lface[6][4] = {{0,1,2,3},
6294 map<NodeTuple<2>,
int> nodesEdgeMap_subdomain;
6297 register int nodes[2];
6301 nodesEdgeMap_subdomain[nt] = i;
6304 map<NodeTuple<4>,
int> nodesBoundaryMap_subdomain;
6308 register int nodes[4];
6314 assert(nodesBoundaryMap_subdomain.find(ebt) == nodesBoundaryMap_subdomain.end());
6315 nodesBoundaryMap_subdomain[ebt] = i;
6320 int local_offset = 0;
6326 if (nN_global_subdomain < nNodes_owned)
6328 subdomain_l2g[eN*nDOF_element + local_offset + nN] = nN_global_subdomain;
6332 subdomain_l2g[eN*nDOF_element + local_offset + nN] = ghostNodeOffset + nN_global_subdomain - nNodes_owned;
6342 for (
int eI=0; eI < 3; eI++)
6343 lagrangeNodesArray[subdomain_l2g[eN*nDOF_element + local_offset +nN]*3+eI] = mesh.
subdomainp->
nodeArray[nN_global_subdomain*3+eI];
6346 for (
int nN = 0; nN < nEdges_element; nN++)
6348 register int nodes[2];
6352 int edge_subdomain = nodesEdgeMap_subdomain[nt];
6353 int edge_global = edgeNumbering_subdomain2global[edge_subdomain];
6354 assert(edge_subdomain >= 0 && edge_global >= 0);
6356 if (edge_subdomain < nEdges_owned)
6357 subdomain_l2g[eN*nDOF_element + local_offset + nN] = nNodes_owned + edge_subdomain;
6359 subdomain_l2g[eN*nDOF_element + local_offset + nN] = ghostEdgeOffset + edge_subdomain - nEdges_owned;
6368 for (
int eI=0; eI < 3; eI++)
6369 lagrangeNodesArray[subdomain_l2g[eN*nDOF_element+local_offset+nN]*3+eI] = 0.5*(mesh.
subdomainp->
nodeArray[nodes[0]*3+eI]+
6373 local_offset += nEdges_element;
6377 register int nodes[4];
6385 int boundary_subdomain = nodesBoundaryMap_subdomain[ebt];
6386 int boundary_global = elementBoundaryNumbering_subdomain2global[boundary_subdomain];
6388 assert(boundary_subdomain >= 0 && boundary_global >= 0);
6391 if (boundary_subdomain < nBoundaries_owned)
6392 subdomain_l2g[eN*nDOF_element + local_offset + nN] = nNodes_owned + nEdges_owned + boundary_subdomain;
6394 subdomain_l2g[eN*nDOF_element + local_offset + nN] = ghostBoundaryOffset + boundary_subdomain - nBoundaries_owned;
6403 for (
int eI=0; eI < 3; eI++)
6404 lagrangeNodesArray[subdomain_l2g[eN*nDOF_element+local_offset+nN]*3+eI] = 0.25*(mesh.
subdomainp->
nodeArray[nodes[0]*3+eI]+
6414 if (eN < nElements_owned)
6415 subdomain_l2g[eN*nDOF_element + local_offset] = nNodes_owned + nEdges_owned + nBoundaries_owned + eN;
6417 subdomain_l2g[eN*nDOF_element + local_offset] = ghostElementOffset + eN - nElements_owned;
6426 for (
int eI=0; eI < 3; eI++)
6439 ISRestoreIndices(quadraticNumberingIS_global_new2old,&quadraticNumbering_global_new2old);
6440 ISDestroy(&quadraticNumberingIS_new2old);
6441 ISDestroy(&quadraticNumberingIS_global_new2old);
6451 const int *elementOffsets_subdomain_owned,
6452 const int *elementNumbering_subdomain2global,
6454 int& nDOF_all_processes,
6455 int& nDOF_subdomain,
6456 int& max_dof_neighbors,
6457 int *offsets_subdomain_owned,
6458 int * subdomain_l2g,
6459 int* subdomain2global)
6462 using namespace std;
6464 ierr = MPI_Comm_size(PROTEUS_COMM_WORLD,&size);
6465 ierr = MPI_Comm_rank(PROTEUS_COMM_WORLD,&rank);
6470 assert(elementOffsets_subdomain_owned);
6471 assert(elementNumbering_subdomain2global);
6474 const int nElements_owned = elementOffsets_subdomain_owned[rank+1]-elementOffsets_subdomain_owned[rank];
6476 assert(offsets_subdomain_owned);
6477 assert(subdomain2global);
6478 assert(subdomain_l2g);
6483 for (
int sdN=0; sdN < size+1; sdN++)
6484 offsets_subdomain_owned[sdN] = elementOffsets_subdomain_owned[sdN]*nDOF_element;
6489 for (
int i = 0; i < nDOF_element; i++)
6491 subdomain2global[eN*nDOF_element + i] = elementNumbering_subdomain2global[eN]*nDOF_element + i;
6494 subdomain_l2g[eN*nDOF_element+i] = eN*nDOF_element + i;