50 #ifndef _ZOLTAN2_PAMGENMESHADAPTER_HPP_ 51 #define _ZOLTAN2_PAMGENMESHADAPTER_HPP_ 55 #include <Teuchos_as.hpp> 59 #include <pamgen_im_exodusII.h> 60 #include <pamgen_im_ne_nemesisI.h> 90 template <
typename User>
136 (
MESH_FACE == etype && 2 == dimension_)) {
150 (
MESH_FACE == etype && 2 == dimension_)) {
151 Ids = element_num_map_;
164 (
MESH_FACE == etype && 2 == dimension_)) {
165 Types = elemTopology;
169 Types = nodeTopology;
176 int &stride,
int idx = 0)
const 185 int &stride,
int dim)
const {
187 (
MESH_FACE == etype && 2 == dimension_)) {
190 }
else if (dim == 1) {
191 coords = Acoords_ + num_elem_;
192 }
else if (dim == 2) {
193 coords = Acoords_ + 2 * num_elem_;
196 }
else if (
MESH_REGION == etype && 2 == dimension_) {
202 }
else if (dim == 1) {
203 coords = coords_ + num_nodes_;
204 }
else if (dim == 2) {
205 coords = coords_ + 2 * num_nodes_;
242 const lno_t *&offsets,
const gno_t *& adjacencyIds)
const 246 offsets = elemOffsets_;
247 adjacencyIds = elemToNode_;
250 offsets = nodeOffsets_;
251 adjacencyIds = nodeToElem_;
252 }
else if (
MESH_REGION == source && 2 == dimension_) {
263 #ifndef USE_MESH_ADAPTER 267 if (sourcetarget ==
MESH_REGION && dimension_ == 3)
return true;
268 if (sourcetarget ==
MESH_FACE && dimension_ == 2)
return true;
271 if (through ==
MESH_REGION && dimension_ == 3)
return true;
272 if (through ==
MESH_FACE && dimension_ == 2)
return true;
281 ((sourcetarget ==
MESH_REGION && dimension_ == 3) ||
282 (sourcetarget ==
MESH_FACE && dimension_ == 2))) {
288 (through ==
MESH_FACE && dimension_ == 2))) {
296 const lno_t *&offsets,
const gno_t *&adjacencyIds)
const 299 ((sourcetarget ==
MESH_REGION && dimension_ == 3) ||
300 (sourcetarget ==
MESH_FACE && dimension_ == 2))) {
302 adjacencyIds = eAdj_;
305 (through ==
MESH_FACE && dimension_ == 2))) {
307 adjacencyIds = nAdj_;
319 (
MESH_FACE == etype && 2 == dimension_) ||
321 return entityDegreeWeight_[idx];
328 int dimension_, num_nodes_global_, num_elems_global_, num_nodes_, num_elem_;
329 gno_t *element_num_map_, *node_num_map_;
331 lno_t tnoct_, *elemOffsets_;
333 lno_t telct_, *nodeOffsets_;
335 int nWeightsPerEntity_;
336 bool *entityDegreeWeight_;
338 scalar_t *coords_, *Acoords_;
339 lno_t *eStart_, *nStart_;
340 gno_t *eAdj_, *nAdj_;
341 size_t nEadj_, nNadj_;
351 template <
typename User>
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 coords_ =
new scalar_t [num_nodes_ * dimension_];
385 error += im_ex_get_coord(exoid, coords_, coords_ + num_nodes_,
386 coords_ + 2 * num_nodes_);
388 element_num_map_ =
new gno_t[num_elem_];
389 std::vector<int> tmp;
390 tmp.resize(num_elem_);
395 error += im_ex_get_elem_num_map(exoid, &tmp[0]);
396 for(
size_t i = 0; i < tmp.size(); i++)
397 element_num_map_[i] = static_cast<gno_t>(tmp[i]);
400 tmp.resize(num_nodes_);
401 node_num_map_ =
new gno_t [num_nodes_];
406 error += im_ex_get_node_num_map(exoid, &tmp[0]);
407 for(
size_t i = 0; i < tmp.size(); i++)
408 node_num_map_[i] = static_cast<gno_t>(tmp[i]);
411 for (
int i=0;i<num_nodes_;i++)
412 nodeTopology[i] =
POINT;
414 for (
int i=0;i<num_elem_;i++) {
421 int *elem_blk_ids =
new int [num_elem_blk];
422 error += im_ex_get_elem_blk_ids(exoid, elem_blk_ids);
424 int *num_nodes_per_elem =
new int [num_elem_blk];
425 int *num_attr =
new int [num_elem_blk];
426 int *num_elem_this_blk =
new int [num_elem_blk];
427 char **elem_type =
new char * [num_elem_blk];
428 int **connect =
new int * [num_elem_blk];
430 for(
int i = 0; i < num_elem_blk; i++){
431 elem_type[i] =
new char [MAX_STR_LENGTH + 1];
432 error += im_ex_get_elem_block(exoid, elem_blk_ids[i], elem_type[i],
433 (
int*)&(num_elem_this_blk[i]),
434 (
int*)&(num_nodes_per_elem[i]),
435 (
int*)&(num_attr[i]));
436 delete[] elem_type[i];
443 Acoords_ =
new scalar_t [num_elem_ * dimension_];
445 std::vector<std::vector<gno_t> > sur_elem;
446 sur_elem.resize(num_nodes_);
448 for(
int b = 0; b < num_elem_blk; b++) {
449 connect[b] =
new int [num_nodes_per_elem[b]*num_elem_this_blk[b]];
450 error += im_ex_get_elem_conn(exoid, elem_blk_ids[b], connect[b]);
452 for(
int i = 0; i < num_elem_this_blk[b]; i++) {
454 Acoords_[num_elem_ + a] = 0;
456 if (3 == dimension_) {
457 Acoords_[2 * num_elem_ + a] = 0;
460 for(
int j = 0; j < num_nodes_per_elem[b]; j++) {
461 int node = connect[b][i * num_nodes_per_elem[b] + j] - 1;
462 Acoords_[a] += coords_[node];
463 Acoords_[num_elem_ + a] += coords_[num_nodes_ + node];
465 if(3 == dimension_) {
466 Acoords_[2 * num_elem_ + a] += coords_[2 * num_nodes_ + node];
475 if (sur_elem[node].empty() ||
476 element_num_map_[a] != sur_elem[node][sur_elem[node].size()-1]) {
478 sur_elem[node].push_back(element_num_map_[a]);
482 Acoords_[a] /= num_nodes_per_elem[b];
483 Acoords_[num_elem_ + a] /= num_nodes_per_elem[b];
485 if(3 == dimension_) {
486 Acoords_[2 * num_elem_ + a] /= num_nodes_per_elem[b];
494 delete[] elem_blk_ids;
496 int nnodes_per_elem = num_nodes_per_elem[0];
497 elemToNode_ =
new gno_t[num_elem_ * nnodes_per_elem];
499 elemOffsets_ =
new lno_t [num_elem_+1];
501 int **reconnect =
new int * [num_elem_];
504 for (
int b = 0; b < num_elem_blk; b++) {
505 for (
int i = 0; i < num_elem_this_blk[b]; i++) {
506 elemOffsets_[telct] = tnoct_;
507 reconnect[telct] =
new int [num_nodes_per_elem[b]];
509 for (
int j = 0; j < num_nodes_per_elem[b]; j++) {
511 as<gno_t>(node_num_map_[connect[b][i*num_nodes_per_elem[b] + j]-1]);
512 reconnect[telct][j] = connect[b][i*num_nodes_per_elem[b] + j];
519 elemOffsets_[telct] = tnoct_;
521 delete[] num_nodes_per_elem;
522 num_nodes_per_elem = NULL;
523 delete[] num_elem_this_blk;
524 num_elem_this_blk = NULL;
526 for(
int b = 0; b < num_elem_blk; b++) {
533 int max_side_nodes = nnodes_per_elem;
534 int *side_nodes =
new int [max_side_nodes];
535 int *mirror_nodes =
new int [max_side_nodes];
538 eStart_ =
new lno_t [num_elem_+1];
539 nStart_ =
new lno_t [num_nodes_+1];
540 std::vector<int> eAdj;
541 std::vector<int> nAdj;
543 for (
int i=0; i < max_side_nodes; i++) {
545 mirror_nodes[i]=-999;
551 for(
int ncnt=0; ncnt < num_nodes_; ncnt++) {
552 if(sur_elem[ncnt].empty()) {
553 std::cout <<
"WARNING: Node = " << ncnt+1 <<
" has no elements" 556 size_t nsur = sur_elem[ncnt].size();
562 nodeToElem_ =
new gno_t[num_nodes_ * max_nsur];
563 nodeOffsets_ =
new lno_t[num_nodes_+1];
566 for (
int ncnt = 0; ncnt < num_nodes_; ncnt++) {
567 nodeOffsets_[ncnt] = telct_;
568 nStart_[ncnt] = nNadj_;
571 for (
size_t i = 0; i < sur_elem[ncnt].size(); i++) {
572 nodeToElem_[telct_] = sur_elem[ncnt][i];
575 #ifndef USE_MESH_ADAPTER 576 for(
int ecnt = 0; ecnt < num_elem_; ecnt++) {
577 if (element_num_map_[ecnt] == sur_elem[ncnt][i]) {
578 for (
int j = 0; j < nnodes_per_elem; j++) {
579 typename MapType::iterator iter =
580 nAdjMap.find(elemToNode_[elemOffsets_[ecnt]+j]);
582 if (node_num_map_[ncnt] != elemToNode_[elemOffsets_[ecnt]+j] &&
583 iter == nAdjMap.end() ) {
584 nAdj.push_back(elemToNode_[elemOffsets_[ecnt]+j]);
586 nAdjMap.insert({elemToNode_[elemOffsets_[ecnt]+j],
587 elemToNode_[elemOffsets_[ecnt]+j]});
600 nodeOffsets_[num_nodes_] = telct_;
601 nStart_[num_nodes_] = nNadj_;
603 nAdj_ =
new gno_t [nNadj_];
605 for (
size_t i=0; i < nNadj_; i++) {
606 nAdj_[i] = as<gno_t>(nAdj[i]);
609 int nprocs = comm.getSize();
611 int neid=0,num_elem_blks_global,num_node_sets_global,num_side_sets_global;
612 error += im_ne_get_init_global(neid,&num_nodes_global_,&num_elems_global_,
613 &num_elem_blks_global,&num_node_sets_global,
614 &num_side_sets_global);
616 int num_internal_nodes, num_border_nodes, num_external_nodes;
617 int num_internal_elems, num_border_elems, num_node_cmaps, num_elem_cmaps;
619 error += im_ne_get_loadbal_param(neid, &num_internal_nodes,
620 &num_border_nodes, &num_external_nodes,
621 &num_internal_elems, &num_border_elems,
622 &num_node_cmaps, &num_elem_cmaps, proc);
624 int *node_cmap_ids =
new int [num_node_cmaps];
625 int *node_cmap_node_cnts =
new int [num_node_cmaps];
626 int *elem_cmap_ids =
new int [num_elem_cmaps];
627 int *elem_cmap_elem_cnts =
new int [num_elem_cmaps];
628 error += im_ne_get_cmap_params(neid, node_cmap_ids, node_cmap_node_cnts,
629 elem_cmap_ids, elem_cmap_elem_cnts, proc);
630 delete[] elem_cmap_ids;
631 elem_cmap_ids = NULL;
632 delete[] elem_cmap_elem_cnts;
633 elem_cmap_elem_cnts = NULL;
635 int **node_ids =
new int * [num_node_cmaps];
636 int **node_proc_ids =
new int * [num_node_cmaps];
637 for(
int j = 0; j < num_node_cmaps; j++) {
638 node_ids[j] =
new int [node_cmap_node_cnts[j]];
639 node_proc_ids[j] =
new int [node_cmap_node_cnts[j]];
640 error += im_ne_get_node_cmap(neid, node_cmap_ids[j], node_ids[j],
641 node_proc_ids[j], proc);
643 delete[] node_cmap_ids;
644 node_cmap_ids = NULL;
645 int *sendCount =
new int [nprocs];
646 int *recvCount =
new int [nprocs];
649 RCP<CommRequest<int> >*requests=
new RCP<CommRequest<int> >[num_node_cmaps];
650 for (
int cnt = 0, i = 0; i < num_node_cmaps; i++) {
653 Teuchos::ireceive<int,int>(comm,
654 rcp(&(recvCount[node_proc_ids[i][0]]),
656 node_proc_ids[i][0]);
661 Teuchos::barrier<int>(comm);
662 size_t totalsend = 0;
664 for(
int j = 0; j < num_node_cmaps; j++) {
665 sendCount[node_proc_ids[j][0]] = 1;
666 for(
int i = 0; i < node_cmap_node_cnts[j]; i++) {
667 sendCount[node_proc_ids[j][i]] += sur_elem[node_ids[j][i]-1].size()+2;
669 totalsend += sendCount[node_proc_ids[j][0]];
673 for (
int i = 0; i < num_node_cmaps; i++) {
675 Teuchos::readySend<int,int>(comm, sendCount[node_proc_ids[i][0]],
676 node_proc_ids[i][0]);
683 Teuchos::waitAll<int>(comm, arrayView(requests, num_node_cmaps));
690 size_t totalrecv = 0;
693 for(
int i = 0; i < num_node_cmaps; i++) {
694 if (recvCount[node_proc_ids[i][0]] > 0) {
695 totalrecv += recvCount[node_proc_ids[i][0]];
697 if (recvCount[node_proc_ids[i][0]] > maxMsg)
698 maxMsg = recvCount[node_proc_ids[i][0]];
703 if (totalrecv) rbuf =
new gno_t[totalrecv];
705 requests =
new RCP<CommRequest<int> > [nrecvranks];
714 if (
size_t(maxMsg) *
sizeof(int) > INT_MAX) OK[0] =
false;
715 if (totalrecv && !rbuf) OK[1] = 0;
716 if (!requests) OK[1] = 0;
722 if (OK[0] && OK[1]) {
724 for (
int i = 0; i < num_node_cmaps; i++) {
725 if (recvCount[node_proc_ids[i][0]]) {
729 ireceive<int,gno_t>(comm,
730 Teuchos::arcp(&rbuf[offset], 0,
731 recvCount[node_proc_ids[i][0]],
733 node_proc_ids[i][0]);
737 offset += recvCount[node_proc_ids[i][0]];
744 Teuchos::reduceAll<int>(comm, Teuchos::REDUCE_MIN, 2, OK, gOK);
745 if (!gOK[0] || !gOK[1]) {
749 throw std::runtime_error(
"Max single message length exceeded");
751 throw std::bad_alloc();
755 if (totalsend) sbuf =
new gno_t[totalsend];
758 for(
int j = 0; j < num_node_cmaps; j++) {
759 sbuf[a++] = node_cmap_node_cnts[j];
760 for(
int i = 0; i < node_cmap_node_cnts[j]; i++) {
761 sbuf[a++] = node_num_map_[node_ids[j][i]-1];
762 sbuf[a++] = sur_elem[node_ids[j][i]-1].size();
763 for(
size_t ecnt=0; ecnt < sur_elem[node_ids[j][i]-1].size(); ecnt++) {
764 sbuf[a++] = sur_elem[node_ids[j][i]-1][ecnt];
769 delete[] node_cmap_node_cnts;
770 node_cmap_node_cnts = NULL;
772 for(
int j = 0; j < num_node_cmaps; j++) {
773 delete[] node_ids[j];
778 ArrayRCP<gno_t> sendBuf;
781 sendBuf = ArrayRCP<gno_t>(sbuf, 0, totalsend,
true);
783 sendBuf = Teuchos::null;
787 for (
int i = 0; i < num_node_cmaps; i++) {
788 if (sendCount[node_proc_ids[i][0]]) {
790 Teuchos::readySend<int, gno_t>(comm,
791 Teuchos::arrayView(&sendBuf[offset],
792 sendCount[node_proc_ids[i][0]]),
793 node_proc_ids[i][0]);
797 offset += sendCount[node_proc_ids[i][0]];
800 for(
int j = 0; j < num_node_cmaps; j++) {
801 delete[] node_proc_ids[j];
804 delete[] node_proc_ids;
805 node_proc_ids = NULL;
810 Teuchos::waitAll<int>(comm, Teuchos::arrayView(requests, nrecvranks));
817 for (
int i = 0; i < num_node_cmaps; i++) {
818 int num_nodes_this_processor = rbuf[a++];
820 for (
int j = 0; j < num_nodes_this_processor; j++) {
821 int this_node = rbuf[a++];
822 int num_elem_this_node = rbuf[a++];
824 for (
int ncnt = 0; ncnt < num_nodes_; ncnt++) {
825 if (node_num_map_[ncnt] == this_node) {
826 for (
int ecnt = 0; ecnt < num_elem_this_node; ecnt++) {
827 sur_elem[ncnt].push_back(rbuf[a++]);
839 #ifndef USE_MESH_ADAPTER 840 for(
int ecnt=0; ecnt < num_elem_; ecnt++) {
841 eStart_[ecnt] = nEadj_;
843 int nnodes = nnodes_per_elem;
844 for(
int ncnt=0; ncnt < nnodes; ncnt++) {
845 int node = reconnect[ecnt][ncnt]-1;
846 for(
size_t i=0; i < sur_elem[node].size(); i++) {
847 int entry = sur_elem[node][i];
848 typename MapType::iterator iter = eAdjMap.find(entry);
850 if(element_num_map_[ecnt] != entry &&
851 iter == eAdjMap.end()) {
852 eAdj.push_back(entry);
854 eAdjMap.insert({entry, entry});
863 for(
int b = 0; b < num_elem_; b++) {
864 delete[] reconnect[b];
869 eStart_[num_elem_] = nEadj_;
871 eAdj_ =
new gno_t [nEadj_];
873 for (
size_t i=0; i < nEadj_; i++) {
874 eAdj_[i] = as<gno_t>(eAdj[i]);
879 delete[] mirror_nodes;
882 if (nWeightsPerEntity_ > 0) {
883 entityDegreeWeight_ =
new bool [nWeightsPerEntity_];
884 for (
int i=0; i < nWeightsPerEntity_; i++) {
885 entityDegreeWeight_[i] =
false;
891 template <
typename User>
894 if (idx >= 0 && idx < nWeightsPerEntity_)
895 entityDegreeWeight_[idx] =
true;
897 std::cout <<
"WARNING: invalid entity weight index, " << idx <<
", ignored" 901 template <
typename User>
904 std::string fn(
" PamgenMesh ");
905 std::cout << me << fn
906 <<
" dim = " << dimension_
907 <<
" nodes = " << num_nodes_
908 <<
" nelems = " << num_elem_
911 for (
int i = 0; i < num_elem_; i++) {
912 std::cout << me << fn << i
913 <<
" Elem " << element_num_map_[i]
915 for (
int j = 0; j < dimension_; j++)
916 std::cout << Acoords_[i + j * num_elem_] <<
" ";
917 std::cout << std::endl;
920 #ifndef USE_MESH_ADAPTER 921 for (
int i = 0; i < num_elem_; i++) {
922 std::cout << me << fn << i
923 <<
" Elem " << element_num_map_[i]
925 for (
int j = eStart_[i]; j < eStart_[i+1]; j++)
926 std::cout << eAdj_[j] <<
" ";
927 std::cout << std::endl;
bool avail2ndAdjs(MeshEntityType sourcetarget, MeshEntityType through) const
Returns whether a second adjacency combination is available. If combination is not available in the M...
InputTraits< User >::part_t part_t
#define Z2_FORWARD_EXCEPTIONS
Forward an exception back through call stack.
void setEntityTypes(std::string ptypestr, std::string atypestr, std::string satypestr)
Sets the primary, adjacency, and second adjacency entity types. Called by algorithm based on paramete...
bool useDegreeAsWeightOf(MeshEntityType etype, int idx) const
Optional method allowing the idx-th weight of entity type etype to be set as the number of neighbors ...
int getDimension() const
Return dimension of the entity coordinates, if any.
Defines the MeshAdapter interface.
MeshAdapter defines the interface for mesh input.
void getWeightsViewOf(MeshEntityType etype, const scalar_t *&weights, int &stride, int idx=0) const
Provide a pointer to one of the number of this process' optional entity weights.
void getIDsViewOf(MeshEntityType etype, const gno_t *&Ids) const
Provide a pointer to this process' identifiers.
InputTraits< User >::gno_t gno_t
bool areEntityIDsUnique(MeshEntityType etype) const
Provide a pointer to the entity topology types.
size_t getLocalNumAdjs(MeshEntityType source, MeshEntityType target) const
Returns the number of first adjacencies on this process.
InputTraits< User >::scalar_t scalar_t
MeshAdapter< User > base_adapter_t
size_t getLocalNum2ndAdjs(MeshEntityType sourcetarget, MeshEntityType through) const
if avail2ndAdjs(), returns the number of second adjacencies on this process.
std::map< gno_t, gno_t > MapType
EntityTopologyType
Enumerate entity topology types for meshes: points,lines,polygons,triangles,quadrilaterals, polyhedrons, tetrahedrons, hexhedrons, prisms, or pyramids.
InputTraits< User >::node_t node_t
void getTopologyViewOf(MeshEntityType etype, enum EntityTopologyType const *&Types) const
Provide a pointer to the entity topology types.
void getCoordinatesViewOf(MeshEntityType etype, const scalar_t *&coords, int &stride, int dim) const
Provide a pointer to one dimension of entity coordinates.
#define Z2_THROW_NOT_IMPLEMENTED_IN_ADAPTER
InputTraits< User >::lno_t lno_t
MeshEntityType
Enumerate entity types for meshes: Regions, Faces, Edges, or Vertices.
void get2ndAdjsView(MeshEntityType sourcetarget, MeshEntityType through, const lno_t *&offsets, const gno_t *&adjacencyIds) const
void getAdjsView(MeshEntityType source, MeshEntityType target, const lno_t *&offsets, const gno_t *&adjacencyIds) const
PamgenMeshAdapter(const Comm< int > &comm, std::string typestr="region", int nEntWgts=0)
Constructor for mesh with identifiers but no coordinates or edges.
void setWeightIsDegree(int idx)
Specify an index for which the weight should be the degree of the entity idx Zoltan2 will use the en...
This class represents a mesh.
size_t getLocalNumOf(MeshEntityType etype) const
Returns the global number of mesh entities of MeshEntityType.
bool availAdjs(MeshEntityType source, MeshEntityType target) const
Returns whether a first adjacency combination is available.
This file defines the StridedData class.