353 std::string typestr,
int nEntWgts):
354 dimension_(0), nWeightsPerEntity_(nEntWgts), entityDegreeWeight_()
361 int num_elem_blk, num_node_sets, num_side_sets;
362 error += im_ex_get_init(exoid, title, &dimension_,
363 &num_nodes_, &num_elem_, &num_elem_blk,
364 &num_node_sets, &num_side_sets);
366 if (typestr.compare(
"region") == 0) {
373 else if (typestr.compare(
"vertex") == 0) {
383 double *dcoords =
new double [num_nodes_ * dimension_];
384 coords_ =
new scalar_t [num_nodes_ * dimension_];
386 error += im_ex_get_coord(exoid, dcoords, dcoords + num_nodes_,
387 dcoords + 2 * num_nodes_);
389 size_t dlen = num_nodes_ * dimension_;
390 for (
size_t i = 0; i < dlen; i++) coords_[i] = as<scalar_t>(dcoords[i]);
393 element_num_map_ =
new gno_t[num_elem_];
394 std::vector<int> tmp;
395 tmp.resize(num_elem_);
400 error += im_ex_get_elem_num_map(exoid, &tmp[0]);
401 for(
size_t i = 0; i < tmp.size(); i++)
402 element_num_map_[i] =
static_cast<gno_t>(tmp[i]);
405 tmp.resize(num_nodes_);
406 node_num_map_ =
new gno_t [num_nodes_];
411 error += im_ex_get_node_num_map(exoid, &tmp[0]);
412 for(
size_t i = 0; i < tmp.size(); i++)
413 node_num_map_[i] =
static_cast<gno_t>(tmp[i]);
416 for (
int i=0;i<num_nodes_;i++)
417 nodeTopology[i] =
POINT;
419 for (
int i=0;i<num_elem_;i++) {
426 int *elem_blk_ids =
new int [num_elem_blk];
427 error += im_ex_get_elem_blk_ids(exoid, elem_blk_ids);
429 int *num_nodes_per_elem =
new int [num_elem_blk];
430 int *num_attr =
new int [num_elem_blk];
431 int *num_elem_this_blk =
new int [num_elem_blk];
432 char **elem_type =
new char * [num_elem_blk];
433 int **connect =
new int * [num_elem_blk];
435 for(
int i = 0; i < num_elem_blk; i++){
436 elem_type[i] =
new char [MAX_STR_LENGTH + 1];
437 error += im_ex_get_elem_block(exoid, elem_blk_ids[i], elem_type[i],
438 (
int*)&(num_elem_this_blk[i]),
439 (
int*)&(num_nodes_per_elem[i]),
440 (
int*)&(num_attr[i]));
441 delete[] elem_type[i];
448 Acoords_ =
new scalar_t [num_elem_ * dimension_];
450 std::vector<std::vector<gno_t> > sur_elem;
451 sur_elem.resize(num_nodes_);
453 for(
int b = 0; b < num_elem_blk; b++) {
454 connect[b] =
new int [num_nodes_per_elem[b]*num_elem_this_blk[b]];
455 error += im_ex_get_elem_conn(exoid, elem_blk_ids[b], connect[b]);
457 for(
int i = 0; i < num_elem_this_blk[b]; i++) {
459 Acoords_[num_elem_ + a] = 0;
461 if (3 == dimension_) {
462 Acoords_[2 * num_elem_ + a] = 0;
465 for(
int j = 0; j < num_nodes_per_elem[b]; j++) {
466 int node = connect[b][i * num_nodes_per_elem[b] + j] - 1;
467 Acoords_[a] += coords_[node];
468 Acoords_[num_elem_ + a] += coords_[num_nodes_ + node];
470 if(3 == dimension_) {
471 Acoords_[2 * num_elem_ + a] += coords_[2 * num_nodes_ + node];
480 if (sur_elem[node].empty() ||
481 element_num_map_[a] != sur_elem[node][sur_elem[node].size()-1]) {
483 sur_elem[node].push_back(element_num_map_[a]);
487 Acoords_[a] /= num_nodes_per_elem[b];
488 Acoords_[num_elem_ + a] /= num_nodes_per_elem[b];
490 if(3 == dimension_) {
491 Acoords_[2 * num_elem_ + a] /= num_nodes_per_elem[b];
499 delete[] elem_blk_ids;
501 int nnodes_per_elem = num_nodes_per_elem[0];
502 elemToNode_ =
new gno_t[num_elem_ * nnodes_per_elem];
504 elemOffsets_ =
new offset_t [num_elem_+1];
506 int **reconnect =
new int * [num_elem_];
509 for (
int b = 0; b < num_elem_blk; b++) {
510 for (
int i = 0; i < num_elem_this_blk[b]; i++) {
511 elemOffsets_[telct] = tnoct_;
512 reconnect[telct] =
new int [num_nodes_per_elem[b]];
514 for (
int j = 0; j < num_nodes_per_elem[b]; j++) {
516 as<gno_t>(node_num_map_[connect[b][i*num_nodes_per_elem[b] + j]-1]);
517 reconnect[telct][j] = connect[b][i*num_nodes_per_elem[b] + j];
524 elemOffsets_[telct] = tnoct_;
526 delete[] num_nodes_per_elem;
527 num_nodes_per_elem = NULL;
528 delete[] num_elem_this_blk;
529 num_elem_this_blk = NULL;
531 for(
int b = 0; b < num_elem_blk; b++) {
538 int max_side_nodes = nnodes_per_elem;
539 int *side_nodes =
new int [max_side_nodes];
540 int *mirror_nodes =
new int [max_side_nodes];
543 eStart_ =
new lno_t [num_elem_+1];
544 nStart_ =
new lno_t [num_nodes_+1];
545 std::vector<int> eAdj;
546 std::vector<int> nAdj;
548 for (
int i=0; i < max_side_nodes; i++) {
550 mirror_nodes[i]=-999;
556 for(
int ncnt=0; ncnt < num_nodes_; ncnt++) {
557 if(sur_elem[ncnt].empty()) {
558 std::cout <<
"WARNING: Node = " << ncnt+1 <<
" has no elements"
561 size_t nsur = sur_elem[ncnt].size();
567 nodeToElem_ =
new gno_t[num_nodes_ * max_nsur];
568 nodeOffsets_ =
new offset_t[num_nodes_+1];
571 for (
int ncnt = 0; ncnt < num_nodes_; ncnt++) {
572 nodeOffsets_[ncnt] = telct_;
573 nStart_[ncnt] = nNadj_;
576 for (
size_t i = 0; i < sur_elem[ncnt].size(); i++) {
577 nodeToElem_[telct_] = sur_elem[ncnt][i];
580#ifndef USE_MESH_ADAPTER
581 for(
int ecnt = 0; ecnt < num_elem_; ecnt++) {
582 if (element_num_map_[ecnt] == sur_elem[ncnt][i]) {
583 for (
int j = 0; j < nnodes_per_elem; j++) {
584 typename MapType::iterator iter =
585 nAdjMap.find(elemToNode_[elemOffsets_[ecnt]+j]);
587 if (node_num_map_[ncnt] != elemToNode_[elemOffsets_[ecnt]+j] &&
588 iter == nAdjMap.end() ) {
589 nAdj.push_back(elemToNode_[elemOffsets_[ecnt]+j]);
591 nAdjMap.insert({elemToNode_[elemOffsets_[ecnt]+j],
592 elemToNode_[elemOffsets_[ecnt]+j]});
605 nodeOffsets_[num_nodes_] = telct_;
606 nStart_[num_nodes_] = nNadj_;
608 nAdj_ =
new gno_t [nNadj_];
610 for (
size_t i=0; i < nNadj_; i++) {
611 nAdj_[i] = as<gno_t>(nAdj[i]);
614 int nprocs = comm.getSize();
616 int neid=0,num_elem_blks_global,num_node_sets_global,num_side_sets_global;
617 error += im_ne_get_init_global(neid,&num_nodes_global_,&num_elems_global_,
618 &num_elem_blks_global,&num_node_sets_global,
619 &num_side_sets_global);
621 int num_internal_nodes, num_border_nodes, num_external_nodes;
622 int num_internal_elems, num_border_elems, num_node_cmaps, num_elem_cmaps;
624 error += im_ne_get_loadbal_param(neid, &num_internal_nodes,
625 &num_border_nodes, &num_external_nodes,
626 &num_internal_elems, &num_border_elems,
627 &num_node_cmaps, &num_elem_cmaps, proc);
629 int *node_cmap_ids =
new int [num_node_cmaps];
630 int *node_cmap_node_cnts =
new int [num_node_cmaps];
631 int *elem_cmap_ids =
new int [num_elem_cmaps];
632 int *elem_cmap_elem_cnts =
new int [num_elem_cmaps];
633 error += im_ne_get_cmap_params(neid, node_cmap_ids, node_cmap_node_cnts,
634 elem_cmap_ids, elem_cmap_elem_cnts, proc);
635 delete[] elem_cmap_ids;
636 elem_cmap_ids = NULL;
637 delete[] elem_cmap_elem_cnts;
638 elem_cmap_elem_cnts = NULL;
640 int **node_ids =
new int * [num_node_cmaps];
641 int **node_proc_ids =
new int * [num_node_cmaps];
642 for(
int j = 0; j < num_node_cmaps; j++) {
643 node_ids[j] =
new int [node_cmap_node_cnts[j]];
644 node_proc_ids[j] =
new int [node_cmap_node_cnts[j]];
645 error += im_ne_get_node_cmap(neid, node_cmap_ids[j], node_ids[j],
646 node_proc_ids[j], proc);
648 delete[] node_cmap_ids;
649 node_cmap_ids = NULL;
650 int *sendCount =
new int [nprocs];
651 int *recvCount =
new int [nprocs];
654 RCP<CommRequest<int> >*requests=
new RCP<CommRequest<int> >[num_node_cmaps];
655 for (
int cnt = 0, i = 0; i < num_node_cmaps; i++) {
658 Teuchos::ireceive<int,int>(comm,
659 rcp(&(recvCount[node_proc_ids[i][0]]),
661 node_proc_ids[i][0]);
666 Teuchos::barrier<int>(comm);
667 size_t totalsend = 0;
669 for(
int j = 0; j < num_node_cmaps; j++) {
670 sendCount[node_proc_ids[j][0]] = 1;
671 for(
int i = 0; i < node_cmap_node_cnts[j]; i++) {
672 sendCount[node_proc_ids[j][i]] += sur_elem[node_ids[j][i]-1].size()+2;
674 totalsend += sendCount[node_proc_ids[j][0]];
678 for (
int i = 0; i < num_node_cmaps; i++) {
680 Teuchos::readySend<int,int>(comm, sendCount[node_proc_ids[i][0]],
681 node_proc_ids[i][0]);
688 Teuchos::waitAll<int>(comm, arrayView(requests, num_node_cmaps));
695 size_t totalrecv = 0;
698 for(
int i = 0; i < num_node_cmaps; i++) {
699 if (recvCount[node_proc_ids[i][0]] > 0) {
700 totalrecv += recvCount[node_proc_ids[i][0]];
702 if (recvCount[node_proc_ids[i][0]] > maxMsg)
703 maxMsg = recvCount[node_proc_ids[i][0]];
708 if (totalrecv) rbuf =
new gno_t[totalrecv];
710 requests =
new RCP<CommRequest<int> > [nrecvranks];
719 if (
size_t(maxMsg) *
sizeof(
int) > INT_MAX) OK[0] =
false;
720 if (totalrecv && !rbuf) OK[1] = 0;
721 if (!requests) OK[1] = 0;
727 if (OK[0] && OK[1]) {
729 for (
int i = 0; i < num_node_cmaps; i++) {
730 if (recvCount[node_proc_ids[i][0]]) {
734 ireceive<int,gno_t>(comm,
735 Teuchos::arcp(&rbuf[offset], 0,
736 recvCount[node_proc_ids[i][0]],
738 node_proc_ids[i][0]);
742 offset += recvCount[node_proc_ids[i][0]];
749 Teuchos::reduceAll<int>(comm, Teuchos::REDUCE_MIN, 2, OK, gOK);
750 if (!gOK[0] || !gOK[1]) {
754 throw std::runtime_error(
"Max single message length exceeded");
756 throw std::bad_alloc();
760 if (totalsend) sbuf =
new gno_t[totalsend];
763 for(
int j = 0; j < num_node_cmaps; j++) {
764 sbuf[a++] = node_cmap_node_cnts[j];
765 for(
int i = 0; i < node_cmap_node_cnts[j]; i++) {
766 sbuf[a++] = node_num_map_[node_ids[j][i]-1];
767 sbuf[a++] = sur_elem[node_ids[j][i]-1].size();
768 for(
size_t ecnt=0; ecnt < sur_elem[node_ids[j][i]-1].size(); ecnt++) {
769 sbuf[a++] = sur_elem[node_ids[j][i]-1][ecnt];
774 delete[] node_cmap_node_cnts;
775 node_cmap_node_cnts = NULL;
777 for(
int j = 0; j < num_node_cmaps; j++) {
778 delete[] node_ids[j];
783 ArrayRCP<gno_t> sendBuf;
786 sendBuf = ArrayRCP<gno_t>(sbuf, 0, totalsend,
true);
788 sendBuf = Teuchos::null;
792 for (
int i = 0; i < num_node_cmaps; i++) {
793 if (sendCount[node_proc_ids[i][0]]) {
795 Teuchos::readySend<int, gno_t>(comm,
796 Teuchos::arrayView(&sendBuf[offset],
797 sendCount[node_proc_ids[i][0]]),
798 node_proc_ids[i][0]);
802 offset += sendCount[node_proc_ids[i][0]];
805 for(
int j = 0; j < num_node_cmaps; j++) {
806 delete[] node_proc_ids[j];
809 delete[] node_proc_ids;
810 node_proc_ids = NULL;
815 Teuchos::waitAll<int>(comm, Teuchos::arrayView(requests, nrecvranks));
822 for (
int i = 0; i < num_node_cmaps; i++) {
823 int num_nodes_this_processor = rbuf[a++];
825 for (
int j = 0; j < num_nodes_this_processor; j++) {
826 int this_node = rbuf[a++];
827 int num_elem_this_node = rbuf[a++];
829 for (
int ncnt = 0; ncnt < num_nodes_; ncnt++) {
830 if (node_num_map_[ncnt] == this_node) {
831 for (
int ecnt = 0; ecnt < num_elem_this_node; ecnt++) {
832 sur_elem[ncnt].push_back(rbuf[a++]);
844#ifndef USE_MESH_ADAPTER
845 for(
int ecnt=0; ecnt < num_elem_; ecnt++) {
846 eStart_[ecnt] = nEadj_;
848 int nnodes = nnodes_per_elem;
849 for(
int ncnt=0; ncnt < nnodes; ncnt++) {
850 int node = reconnect[ecnt][ncnt]-1;
851 for(
size_t i=0; i < sur_elem[node].size(); i++) {
852 int entry = sur_elem[node][i];
853 typename MapType::iterator iter = eAdjMap.find(entry);
855 if(element_num_map_[ecnt] != entry &&
856 iter == eAdjMap.end()) {
857 eAdj.push_back(entry);
859 eAdjMap.insert({entry, entry});
868 for(
int b = 0; b < num_elem_; b++) {
869 delete[] reconnect[b];
874 eStart_[num_elem_] = nEadj_;
876 eAdj_ =
new gno_t [nEadj_];
878 for (
size_t i=0; i < nEadj_; i++) {
879 eAdj_[i] = as<gno_t>(eAdj[i]);
884 delete[] mirror_nodes;
887 if (nWeightsPerEntity_ > 0) {
888 entityDegreeWeight_ =
new bool [nWeightsPerEntity_];
889 for (
int i=0; i < nWeightsPerEntity_; i++) {
890 entityDegreeWeight_[i] =
false;