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>
135 (
MESH_FACE == etype && 2 == dimension_)) {
149 (
MESH_FACE == etype && 2 == dimension_)) {
150 Ids = element_num_map_;
163 (
MESH_FACE == etype && 2 == dimension_)) {
164 Types = elemTopology;
168 Types = nodeTopology;
175 int &stride,
int idx = 0)
const 184 int &stride,
int dim)
const {
186 (
MESH_FACE == etype && 2 == dimension_)) {
189 }
else if (dim == 1) {
190 coords = Acoords_ + num_elem_;
191 }
else if (dim == 2) {
192 coords = Acoords_ + 2 * num_elem_;
195 }
else if (
MESH_REGION == etype && 2 == dimension_) {
201 }
else if (dim == 1) {
202 coords = coords_ + num_nodes_;
203 }
else if (dim == 2) {
204 coords = coords_ + 2 * num_nodes_;
241 const lno_t *&offsets,
const gno_t *& adjacencyIds)
const 245 offsets = elemOffsets_;
246 adjacencyIds = elemToNode_;
249 offsets = nodeOffsets_;
250 adjacencyIds = nodeToElem_;
251 }
else if (
MESH_REGION == source && 2 == dimension_) {
262 #ifndef USE_MESH_ADAPTER 266 if (sourcetarget ==
MESH_REGION && dimension_ == 3)
return true;
267 if (sourcetarget ==
MESH_FACE && dimension_ == 2)
return true;
270 if (through ==
MESH_REGION && dimension_ == 3)
return true;
271 if (through ==
MESH_FACE && dimension_ == 2)
return true;
280 ((sourcetarget ==
MESH_REGION && dimension_ == 3) ||
281 (sourcetarget ==
MESH_FACE && dimension_ == 2))) {
287 (through ==
MESH_FACE && dimension_ == 2))) {
295 const lno_t *&offsets,
const gno_t *&adjacencyIds)
const 298 ((sourcetarget ==
MESH_REGION && dimension_ == 3) ||
299 (sourcetarget ==
MESH_FACE && dimension_ == 2))) {
301 adjacencyIds = eAdj_;
304 (through ==
MESH_FACE && dimension_ == 2))) {
306 adjacencyIds = nAdj_;
318 (
MESH_FACE == etype && 2 == dimension_) ||
320 return entityDegreeWeight_[idx];
327 int dimension_, num_nodes_global_, num_elems_global_, num_nodes_, num_elem_;
328 gno_t *element_num_map_, *node_num_map_;
330 lno_t tnoct_, *elemOffsets_;
332 lno_t telct_, *nodeOffsets_;
334 int nWeightsPerEntity_;
335 bool *entityDegreeWeight_;
338 lno_t *eStart_, *nStart_;
339 gno_t *eAdj_, *nAdj_;
340 size_t nEadj_, nNadj_;
350 template <
typename User>
352 std::string typestr,
int nEntWgts):
353 dimension_(0), nWeightsPerEntity_(nEntWgts), entityDegreeWeight_()
360 int num_elem_blk, num_node_sets, num_side_sets;
361 error += im_ex_get_init(exoid, title, &dimension_,
362 &num_nodes_, &num_elem_, &num_elem_blk,
363 &num_node_sets, &num_side_sets);
365 if (typestr.compare(
"region") == 0) {
372 else if (typestr.compare(
"vertex") == 0) {
382 coords_ =
new scalar_t [num_nodes_ * dimension_];
384 error += im_ex_get_coord(exoid, coords_, coords_ + num_nodes_,
385 coords_ + 2 * num_nodes_);
387 element_num_map_ =
new gno_t[num_elem_];
388 std::vector<int> tmp;
389 tmp.resize(num_elem_);
394 error += im_ex_get_elem_num_map(exoid, &tmp[0]);
395 for(
size_t i = 0; i < tmp.size(); i++)
396 element_num_map_[i] = static_cast<gno_t>(tmp[i]);
399 tmp.resize(num_nodes_);
400 node_num_map_ =
new gno_t [num_nodes_];
405 error += im_ex_get_node_num_map(exoid, &tmp[0]);
406 for(
size_t i = 0; i < tmp.size(); i++)
407 node_num_map_[i] = static_cast<gno_t>(tmp[i]);
410 for (
int i=0;i<num_nodes_;i++)
411 nodeTopology[i] =
POINT;
413 for (
int i=0;i<num_elem_;i++) {
420 int *elem_blk_ids =
new int [num_elem_blk];
421 error += im_ex_get_elem_blk_ids(exoid, elem_blk_ids);
423 int *num_nodes_per_elem =
new int [num_elem_blk];
424 int *num_attr =
new int [num_elem_blk];
425 int *num_elem_this_blk =
new int [num_elem_blk];
426 char **elem_type =
new char * [num_elem_blk];
427 int **connect =
new int * [num_elem_blk];
429 for(
int i = 0; i < num_elem_blk; i++){
430 elem_type[i] =
new char [MAX_STR_LENGTH + 1];
431 error += im_ex_get_elem_block(exoid, elem_blk_ids[i], elem_type[i],
432 (
int*)&(num_elem_this_blk[i]),
433 (
int*)&(num_nodes_per_elem[i]),
434 (
int*)&(num_attr[i]));
435 delete[] elem_type[i];
442 Acoords_ =
new scalar_t [num_elem_ * dimension_];
444 std::vector<std::vector<gno_t> > sur_elem;
445 sur_elem.resize(num_nodes_);
447 for(
int b = 0; b < num_elem_blk; b++) {
448 connect[b] =
new int [num_nodes_per_elem[b]*num_elem_this_blk[b]];
449 error += im_ex_get_elem_conn(exoid, elem_blk_ids[b], connect[b]);
451 for(
int i = 0; i < num_elem_this_blk[b]; i++) {
453 Acoords_[num_elem_ + a] = 0;
455 if (3 == dimension_) {
456 Acoords_[2 * num_elem_ + a] = 0;
459 for(
int j = 0; j < num_nodes_per_elem[b]; j++) {
460 int node = connect[b][i * num_nodes_per_elem[b] + j] - 1;
461 Acoords_[a] += coords_[node];
462 Acoords_[num_elem_ + a] += coords_[num_nodes_ + node];
464 if(3 == dimension_) {
465 Acoords_[2 * num_elem_ + a] += coords_[2 * num_nodes_ + node];
474 if (sur_elem[node].empty() ||
475 element_num_map_[a] != sur_elem[node][sur_elem[node].size()-1]) {
477 sur_elem[node].push_back(element_num_map_[a]);
481 Acoords_[a] /= num_nodes_per_elem[b];
482 Acoords_[num_elem_ + a] /= num_nodes_per_elem[b];
484 if(3 == dimension_) {
485 Acoords_[2 * num_elem_ + a] /= num_nodes_per_elem[b];
493 delete[] elem_blk_ids;
495 int nnodes_per_elem = num_nodes_per_elem[0];
496 elemToNode_ =
new gno_t[num_elem_ * nnodes_per_elem];
498 elemOffsets_ =
new lno_t [num_elem_+1];
500 int **reconnect =
new int * [num_elem_];
503 for (
int b = 0; b < num_elem_blk; b++) {
504 for (
int i = 0; i < num_elem_this_blk[b]; i++) {
505 elemOffsets_[telct] = tnoct_;
506 reconnect[telct] =
new int [num_nodes_per_elem[b]];
508 for (
int j = 0; j < num_nodes_per_elem[b]; j++) {
510 as<gno_t>(node_num_map_[connect[b][i*num_nodes_per_elem[b] + j]-1]);
511 reconnect[telct][j] = connect[b][i*num_nodes_per_elem[b] + j];
518 elemOffsets_[telct] = tnoct_;
520 delete[] num_nodes_per_elem;
521 num_nodes_per_elem = NULL;
522 delete[] num_elem_this_blk;
523 num_elem_this_blk = NULL;
525 for(
int b = 0; b < num_elem_blk; b++) {
532 int max_side_nodes = nnodes_per_elem;
533 int *side_nodes =
new int [max_side_nodes];
534 int *mirror_nodes =
new int [max_side_nodes];
537 eStart_ =
new lno_t [num_elem_+1];
538 nStart_ =
new lno_t [num_nodes_+1];
539 std::vector<int> eAdj;
540 std::vector<int> nAdj;
542 for (
int i=0; i < max_side_nodes; i++) {
544 mirror_nodes[i]=-999;
550 for(
int ncnt=0; ncnt < num_nodes_; ncnt++) {
551 if(sur_elem[ncnt].empty()) {
552 std::cout <<
"WARNING: Node = " << ncnt+1 <<
" has no elements" 555 size_t nsur = sur_elem[ncnt].size();
561 nodeToElem_ =
new gno_t[num_nodes_ * max_nsur];
562 nodeOffsets_ =
new lno_t[num_nodes_+1];
565 for (
int ncnt = 0; ncnt < num_nodes_; ncnt++) {
566 nodeOffsets_[ncnt] = telct_;
567 nStart_[ncnt] = nNadj_;
570 for (
size_t i = 0; i < sur_elem[ncnt].size(); i++) {
571 nodeToElem_[telct_] = sur_elem[ncnt][i];
574 #ifndef USE_MESH_ADAPTER 575 for(
int ecnt = 0; ecnt < num_elem_; ecnt++) {
576 if (element_num_map_[ecnt] == sur_elem[ncnt][i]) {
577 for (
int j = 0; j < nnodes_per_elem; j++) {
578 typename MapType::iterator iter =
579 nAdjMap.find(elemToNode_[elemOffsets_[ecnt]+j]);
581 if (node_num_map_[ncnt] != elemToNode_[elemOffsets_[ecnt]+j] &&
582 iter == nAdjMap.end() ) {
583 nAdj.push_back(elemToNode_[elemOffsets_[ecnt]+j]);
585 nAdjMap.insert({elemToNode_[elemOffsets_[ecnt]+j],
586 elemToNode_[elemOffsets_[ecnt]+j]});
599 nodeOffsets_[num_nodes_] = telct_;
600 nStart_[num_nodes_] = nNadj_;
602 nAdj_ =
new gno_t [nNadj_];
604 for (
size_t i=0; i < nNadj_; i++) {
605 nAdj_[i] = as<gno_t>(nAdj[i]);
608 int nprocs = comm.getSize();
610 int neid=0,num_elem_blks_global,num_node_sets_global,num_side_sets_global;
611 error += im_ne_get_init_global(neid,&num_nodes_global_,&num_elems_global_,
612 &num_elem_blks_global,&num_node_sets_global,
613 &num_side_sets_global);
615 int num_internal_nodes, num_border_nodes, num_external_nodes;
616 int num_internal_elems, num_border_elems, num_node_cmaps, num_elem_cmaps;
618 error += im_ne_get_loadbal_param(neid, &num_internal_nodes,
619 &num_border_nodes, &num_external_nodes,
620 &num_internal_elems, &num_border_elems,
621 &num_node_cmaps, &num_elem_cmaps, proc);
623 int *node_cmap_ids =
new int [num_node_cmaps];
624 int *node_cmap_node_cnts =
new int [num_node_cmaps];
625 int *elem_cmap_ids =
new int [num_elem_cmaps];
626 int *elem_cmap_elem_cnts =
new int [num_elem_cmaps];
627 error += im_ne_get_cmap_params(neid, node_cmap_ids, node_cmap_node_cnts,
628 elem_cmap_ids, elem_cmap_elem_cnts, proc);
629 delete[] elem_cmap_ids;
630 elem_cmap_ids = NULL;
631 delete[] elem_cmap_elem_cnts;
632 elem_cmap_elem_cnts = NULL;
634 int **node_ids =
new int * [num_node_cmaps];
635 int **node_proc_ids =
new int * [num_node_cmaps];
636 for(
int j = 0; j < num_node_cmaps; j++) {
637 node_ids[j] =
new int [node_cmap_node_cnts[j]];
638 node_proc_ids[j] =
new int [node_cmap_node_cnts[j]];
639 error += im_ne_get_node_cmap(neid, node_cmap_ids[j], node_ids[j],
640 node_proc_ids[j], proc);
642 delete[] node_cmap_ids;
643 node_cmap_ids = NULL;
644 int *sendCount =
new int [nprocs];
645 int *recvCount =
new int [nprocs];
648 RCP<CommRequest<int> >*requests=
new RCP<CommRequest<int> >[num_node_cmaps];
649 for (
int cnt = 0, i = 0; i < num_node_cmaps; i++) {
652 Teuchos::ireceive<int,int>(comm,
653 rcp(&(recvCount[node_proc_ids[i][0]]),
655 node_proc_ids[i][0]);
660 Teuchos::barrier<int>(comm);
661 size_t totalsend = 0;
663 for(
int j = 0; j < num_node_cmaps; j++) {
664 sendCount[node_proc_ids[j][0]] = 1;
665 for(
int i = 0; i < node_cmap_node_cnts[j]; i++) {
666 sendCount[node_proc_ids[j][i]] += sur_elem[node_ids[j][i]-1].size()+2;
668 totalsend += sendCount[node_proc_ids[j][0]];
672 for (
int i = 0; i < num_node_cmaps; i++) {
674 Teuchos::readySend<int,int>(comm, sendCount[node_proc_ids[i][0]],
675 node_proc_ids[i][0]);
682 Teuchos::waitAll<int>(comm, arrayView(requests, num_node_cmaps));
689 size_t totalrecv = 0;
692 for(
int i = 0; i < num_node_cmaps; i++) {
693 if (recvCount[node_proc_ids[i][0]] > 0) {
694 totalrecv += recvCount[node_proc_ids[i][0]];
696 if (recvCount[node_proc_ids[i][0]] > maxMsg)
697 maxMsg = recvCount[node_proc_ids[i][0]];
702 if (totalrecv) rbuf =
new gno_t[totalrecv];
704 requests =
new RCP<CommRequest<int> > [nrecvranks];
713 if (
size_t(maxMsg) *
sizeof(int) > INT_MAX) OK[0] =
false;
714 if (totalrecv && !rbuf) OK[1] = 0;
715 if (!requests) OK[1] = 0;
721 if (OK[0] && OK[1]) {
723 for (
int i = 0; i < num_node_cmaps; i++) {
724 if (recvCount[node_proc_ids[i][0]]) {
728 ireceive<int,gno_t>(comm,
729 Teuchos::arcp(&rbuf[offset], 0,
730 recvCount[node_proc_ids[i][0]],
732 node_proc_ids[i][0]);
736 offset += recvCount[node_proc_ids[i][0]];
743 Teuchos::reduceAll<int>(comm, Teuchos::REDUCE_MIN, 2, OK, gOK);
744 if (!gOK[0] || !gOK[1]) {
748 throw std::runtime_error(
"Max single message length exceeded");
750 throw std::bad_alloc();
754 if (totalsend) sbuf =
new gno_t[totalsend];
757 for(
int j = 0; j < num_node_cmaps; j++) {
758 sbuf[a++] = node_cmap_node_cnts[j];
759 for(
int i = 0; i < node_cmap_node_cnts[j]; i++) {
760 sbuf[a++] = node_num_map_[node_ids[j][i]-1];
761 sbuf[a++] = sur_elem[node_ids[j][i]-1].size();
762 for(
size_t ecnt=0; ecnt < sur_elem[node_ids[j][i]-1].size(); ecnt++) {
763 sbuf[a++] = sur_elem[node_ids[j][i]-1][ecnt];
768 delete[] node_cmap_node_cnts;
769 node_cmap_node_cnts = NULL;
771 for(
int j = 0; j < num_node_cmaps; j++) {
772 delete[] node_ids[j];
777 ArrayRCP<gno_t> sendBuf;
780 sendBuf = ArrayRCP<gno_t>(sbuf, 0, totalsend,
true);
782 sendBuf = Teuchos::null;
786 for (
int i = 0; i < num_node_cmaps; i++) {
787 if (sendCount[node_proc_ids[i][0]]) {
789 Teuchos::readySend<int, gno_t>(comm,
790 Teuchos::arrayView(&sendBuf[offset],
791 sendCount[node_proc_ids[i][0]]),
792 node_proc_ids[i][0]);
796 offset += sendCount[node_proc_ids[i][0]];
799 for(
int j = 0; j < num_node_cmaps; j++) {
800 delete[] node_proc_ids[j];
803 delete[] node_proc_ids;
804 node_proc_ids = NULL;
809 Teuchos::waitAll<int>(comm, Teuchos::arrayView(requests, nrecvranks));
816 for (
int i = 0; i < num_node_cmaps; i++) {
817 int num_nodes_this_processor = rbuf[a++];
819 for (
int j = 0; j < num_nodes_this_processor; j++) {
820 int this_node = rbuf[a++];
821 int num_elem_this_node = rbuf[a++];
823 for (
int ncnt = 0; ncnt < num_nodes_; ncnt++) {
824 if (node_num_map_[ncnt] == this_node) {
825 for (
int ecnt = 0; ecnt < num_elem_this_node; ecnt++) {
826 sur_elem[ncnt].push_back(rbuf[a++]);
838 #ifndef USE_MESH_ADAPTER 839 for(
int ecnt=0; ecnt < num_elem_; ecnt++) {
840 eStart_[ecnt] = nEadj_;
842 int nnodes = nnodes_per_elem;
843 for(
int ncnt=0; ncnt < nnodes; ncnt++) {
844 int node = reconnect[ecnt][ncnt]-1;
845 for(
size_t i=0; i < sur_elem[node].size(); i++) {
846 int entry = sur_elem[node][i];
847 typename MapType::iterator iter = eAdjMap.find(entry);
849 if(element_num_map_[ecnt] != entry &&
850 iter == eAdjMap.end()) {
851 eAdj.push_back(entry);
853 eAdjMap.insert({entry, entry});
862 for(
int b = 0; b < num_elem_; b++) {
863 delete[] reconnect[b];
868 eStart_[num_elem_] = nEadj_;
870 eAdj_ =
new gno_t [nEadj_];
872 for (
size_t i=0; i < nEadj_; i++) {
873 eAdj_[i] = as<gno_t>(eAdj[i]);
878 delete[] mirror_nodes;
881 if (nWeightsPerEntity_ > 0) {
882 entityDegreeWeight_ =
new bool [nWeightsPerEntity_];
883 for (
int i=0; i < nWeightsPerEntity_; i++) {
884 entityDegreeWeight_[i] =
false;
890 template <
typename User>
893 if (idx >= 0 && idx < nWeightsPerEntity_)
894 entityDegreeWeight_[idx] =
true;
896 std::cout <<
"WARNING: invalid entity weight index, " << idx <<
", ignored" 900 template <
typename User>
903 std::string fn(
" PamgenMesh ");
904 std::cout << me << fn
905 <<
" dim = " << dimension_
906 <<
" nodes = " << num_nodes_
907 <<
" nelems = " << num_elem_
910 for (
int i = 0; i < num_elem_; i++) {
911 std::cout << me << fn << i
912 <<
" Elem " << element_num_map_[i]
914 for (
int j = 0; j < dimension_; j++)
915 std::cout << Acoords_[i + j * num_elem_] <<
" ";
916 std::cout << std::endl;
919 #ifndef USE_MESH_ADAPTER 920 for (
int i = 0; i < num_elem_; i++) {
921 std::cout << me << fn << i
922 <<
" Elem " << element_num_map_[i]
924 for (
int j = eStart_[i]; j < eStart_[i+1]; j++)
925 std::cout << eAdj_[j] <<
" ";
926 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...
size_t getLocalNumOf(MeshEntityType etype) const
Returns the global number of mesh entities of MeshEntityType.
Defines the MeshAdapter interface.
MeshAdapter defines the interface for mesh input.
void getTopologyViewOf(MeshEntityType etype, enum EntityTopologyType const *&Types) const
Provide a pointer to the entity topology types.
#define Z2_THROW_NOT_IMPLEMENTED
InputTraits< User >::gno_t gno_t
InputTraits< User >::scalar_t scalar_t
bool availAdjs(MeshEntityType source, MeshEntityType target) const
Returns whether a first adjacency combination is available.
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 ...
void get2ndAdjsView(MeshEntityType sourcetarget, MeshEntityType through, const lno_t *&offsets, const gno_t *&adjacencyIds) const
if avail2ndAdjs(), set pointers to this process' second adjacencies
void getIDsViewOf(MeshEntityType etype, const gno_t *&Ids) const
Provide a pointer to this process' identifiers.
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 getCoordinatesViewOf(MeshEntityType etype, const scalar_t *&coords, int &stride, int dim) const
Provide a pointer to one dimension of entity coordinates.
size_t getLocalNum2ndAdjs(MeshEntityType sourcetarget, MeshEntityType through) const
if avail2ndAdjs(), returns the number of second adjacencies on this process.
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.
size_t getLocalNumAdjs(MeshEntityType source, MeshEntityType target) const
Returns the number of first adjacencies on this process.
InputTraits< User >::lno_t lno_t
MeshEntityType
Enumerate entity types for meshes: Regions, Faces, Edges, or Vertices.
bool areEntityIDsUnique(MeshEntityType etype) const
Provide a pointer to the entity topology types.
PamgenMeshAdapter(const Comm< int > &comm, std::string typestr="region", int nEntWgts=0)
Constructor for mesh with identifiers but no coordinates or edges.
int getDimension() const
Return dimension of the entity coordinates, if any.
void setWeightIsDegree(int idx)
Specify an index for which the weight should be the degree of the entity idx Zoltan2 will use the en...
void getAdjsView(MeshEntityType source, MeshEntityType target, const lno_t *&offsets, const gno_t *&adjacencyIds) const
Sets pointers to this process' mesh first adjacencies.
This class represents a mesh.
This file defines the StridedData class.