Zoltan2
Zoltan2_PamgenMeshAdapter.hpp
Go to the documentation of this file.
1 // @HEADER
2 //
3 // ***********************************************************************
4 //
5 // Zoltan2: A package of combinatorial algorithms for scientific computing
6 // Copyright 2012 Sandia Corporation
7 //
8 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
9 // the U.S. Government retains certain rights in this software.
10 //
11 // Redistribution and use in source and binary forms, with or without
12 // modification, are permitted provided that the following conditions are
13 // met:
14 //
15 // 1. Redistributions of source code must retain the above copyright
16 // notice, this list of conditions and the following disclaimer.
17 //
18 // 2. Redistributions in binary form must reproduce the above copyright
19 // notice, this list of conditions and the following disclaimer in the
20 // documentation and/or other materials provided with the distribution.
21 //
22 // 3. Neither the name of the Corporation nor the names of the
23 // contributors may be used to endorse or promote products derived from
24 // this software without specific prior written permission.
25 //
26 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
27 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
28 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
29 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
30 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
31 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
32 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
33 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
34 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
35 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
36 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
37 //
38 // Questions? Contact Karen Devine (kddevin@sandia.gov)
39 // Erik Boman (egboman@sandia.gov)
40 // Siva Rajamanickam (srajama@sandia.gov)
41 //
42 // ***********************************************************************
43 //
44 // @HEADER
45 
50 #ifndef _ZOLTAN2_PAMGENMESHADAPTER_HPP_
51 #define _ZOLTAN2_PAMGENMESHADAPTER_HPP_
52 
53 #include <Zoltan2_MeshAdapter.hpp>
54 #include <Zoltan2_StridedData.hpp>
55 #include <Teuchos_as.hpp>
56 #include <vector>
57 #include <string>
58 
59 #include <pamgen_im_exodusII.h>
60 #include <pamgen_im_ne_nemesisI.h>
61 
62 namespace Zoltan2 {
63 
90 template <typename User>
91  class PamgenMeshAdapter: public MeshAdapter<User> {
92 
93 public:
94 
96  typedef typename InputTraits<User>::lno_t lno_t;
97  typedef typename InputTraits<User>::gno_t gno_t;
100  typedef User user_t;
101  typedef std::map<gno_t, gno_t> MapType;
102 
110  PamgenMeshAdapter(const Comm<int> &comm, std::string typestr="region",
111  int nEntWgts=0);
112 
119  void setWeightIsDegree(int idx);
120 
121  void print(int);
122 
124  // The MeshAdapter interface.
125  // This is the interface that would be called by a model or a problem .
127 
128  bool areEntityIDsUnique(MeshEntityType etype) const {
129  return etype==MESH_REGION;
130  }
131 
132  size_t getLocalNumOf(MeshEntityType etype) const
133  {
134  if ((MESH_REGION == etype && 3 == dimension_) ||
135  (MESH_FACE == etype && 2 == dimension_)) {
136  return num_elem_;
137  }
138 
139  if (MESH_VERTEX == etype) {
140  return num_nodes_;
141  }
142 
143  return 0;
144  }
145 
146  void getIDsViewOf(MeshEntityType etype, const gno_t *&Ids) const
147  {
148  if ((MESH_REGION == etype && 3 == dimension_) ||
149  (MESH_FACE == etype && 2 == dimension_)) {
150  Ids = element_num_map_;
151  }
152 
153  else if (MESH_VERTEX == etype) {
154  Ids = node_num_map_;
155  }
156 
157  else Ids = NULL;
158  }
159 
161  enum EntityTopologyType const *&Types) const {
162  if ((MESH_REGION == etype && 3 == dimension_) ||
163  (MESH_FACE == etype && 2 == dimension_)) {
164  Types = elemTopology;
165  }
166 
167  else if (MESH_VERTEX == etype) {
168  Types = nodeTopology;
169  }
170 
171  else Types = NULL;
172  }
173 
175  int &stride, int idx = 0) const
176  {
177  weights = NULL;
178  stride = 0;
179  }
180 
181  int getDimension() const { return dimension_; }
182 
183  void getCoordinatesViewOf(MeshEntityType etype, const scalar_t *&coords,
184  int &stride, int dim) const {
185  if ((MESH_REGION == etype && 3 == dimension_) ||
186  (MESH_FACE == etype && 2 == dimension_)) {
187  if (dim == 0) {
188  coords = Acoords_;
189  } else if (dim == 1) {
190  coords = Acoords_ + num_elem_;
191  } else if (dim == 2) {
192  coords = Acoords_ + 2 * num_elem_;
193  }
194  stride = 1;
195  } else if (MESH_REGION == etype && 2 == dimension_) {
196  coords = NULL;
197  stride = 0;
198  } else if (MESH_VERTEX == etype) {
199  if (dim == 0) {
200  coords = coords_;
201  } else if (dim == 1) {
202  coords = coords_ + num_nodes_;
203  } else if (dim == 2) {
204  coords = coords_ + 2 * num_nodes_;
205  }
206  stride = 1;
207  } else {
208  coords = NULL;
209  stride = 0;
211  }
212  }
213 
214  bool availAdjs(MeshEntityType source, MeshEntityType target) const {
215  if ((MESH_REGION == source && MESH_VERTEX == target && 3 == dimension_) ||
216  (MESH_FACE == source && MESH_VERTEX == target && 2 == dimension_) ||
217  (MESH_VERTEX == source && MESH_REGION == target && 3 == dimension_) ||
218  (MESH_VERTEX == source && MESH_FACE == target && 2 == dimension_)) {
219  return TRUE;
220  }
221 
222  return false;
223  }
224 
225  size_t getLocalNumAdjs(MeshEntityType source, MeshEntityType target) const
226  {
227  if ((MESH_REGION == source && MESH_VERTEX == target && 3 == dimension_) ||
228  (MESH_FACE == source && MESH_VERTEX == target && 2 == dimension_)) {
229  return tnoct_;
230  }
231 
232  if ((MESH_VERTEX == source && MESH_REGION == target && 3 == dimension_) ||
233  (MESH_VERTEX == source && MESH_FACE == target && 2 == dimension_)) {
234  return telct_;
235  }
236 
237  return 0;
238  }
239 
241  const lno_t *&offsets, const gno_t *& adjacencyIds) const
242  {
243  if ((MESH_REGION == source && MESH_VERTEX == target && 3 == dimension_) ||
244  (MESH_FACE == source && MESH_VERTEX == target && 2 == dimension_)) {
245  offsets = elemOffsets_;
246  adjacencyIds = elemToNode_;
247  } else if ((MESH_REGION==target && MESH_VERTEX==source && 3==dimension_) ||
248  (MESH_FACE==target && MESH_VERTEX==source && 2==dimension_)) {
249  offsets = nodeOffsets_;
250  adjacencyIds = nodeToElem_;
251  } else if (MESH_REGION == source && 2 == dimension_) {
252  offsets = NULL;
253  adjacencyIds = NULL;
254  } else {
255  offsets = NULL;
256  adjacencyIds = NULL;
258  }
259  }
260 
261  //#define USE_MESH_ADAPTER
262 #ifndef USE_MESH_ADAPTER
263  bool avail2ndAdjs(MeshEntityType sourcetarget, MeshEntityType through) const
264  {
265  if (through == MESH_VERTEX) {
266  if (sourcetarget == MESH_REGION && dimension_ == 3) return true;
267  if (sourcetarget == MESH_FACE && dimension_ == 2) return true;
268  }
269  if (sourcetarget == MESH_VERTEX) {
270  if (through == MESH_REGION && dimension_ == 3) return true;
271  if (through == MESH_FACE && dimension_ == 2) return true;
272  }
273  return false;
274  }
275 
276  size_t getLocalNum2ndAdjs(MeshEntityType sourcetarget,
277  MeshEntityType through) const
278  {
279  if (through == MESH_VERTEX &&
280  ((sourcetarget == MESH_REGION && dimension_ == 3) ||
281  (sourcetarget == MESH_FACE && dimension_ == 2))) {
282  return nEadj_;
283  }
284 
285  if (sourcetarget == MESH_VERTEX &&
286  ((through == MESH_REGION && dimension_ == 3) ||
287  (through == MESH_FACE && dimension_ == 2))) {
288  return nNadj_;
289  }
290 
291  return 0;
292  }
293 
294  void get2ndAdjsView(MeshEntityType sourcetarget, MeshEntityType through,
295  const lno_t *&offsets, const gno_t *&adjacencyIds) const
296  {
297  if (through == MESH_VERTEX &&
298  ((sourcetarget == MESH_REGION && dimension_ == 3) ||
299  (sourcetarget == MESH_FACE && dimension_ == 2))) {
300  offsets = eStart_;
301  adjacencyIds = eAdj_;
302  } else if (sourcetarget == MESH_VERTEX &&
303  ((through == MESH_REGION && dimension_ == 3) ||
304  (through == MESH_FACE && dimension_ == 2))) {
305  offsets = nStart_;
306  adjacencyIds = nAdj_;
307  } else {
308  offsets = NULL;
309  adjacencyIds = NULL;
311  }
312  }
313 #endif
314 
315  bool useDegreeAsWeightOf(MeshEntityType etype, int idx) const
316  {
317  if ((MESH_REGION == etype && 3 == dimension_) ||
318  (MESH_FACE == etype && 2 == dimension_) ||
319  (MESH_VERTEX == etype)) {
320  return entityDegreeWeight_[idx];
321  }
322 
323  return false;
324  }
325 
326 private:
327  int dimension_, num_nodes_global_, num_elems_global_, num_nodes_, num_elem_;
328  gno_t *element_num_map_, *node_num_map_;
329  gno_t *elemToNode_;
330  lno_t tnoct_, *elemOffsets_;
331  gno_t *nodeToElem_;
332  lno_t telct_, *nodeOffsets_;
333 
334  int nWeightsPerEntity_;
335  bool *entityDegreeWeight_;
336 
337  scalar_t *coords_, *Acoords_;
338  lno_t *eStart_, *nStart_;
339  gno_t *eAdj_, *nAdj_;
340  size_t nEadj_, nNadj_;
341  EntityTopologyType* nodeTopology;
342  EntityTopologyType* elemTopology;
343 
344 };
345 
347 // Definitions
349 
350 template <typename User>
352  std::string typestr, int nEntWgts):
353  dimension_(0), nWeightsPerEntity_(nEntWgts), entityDegreeWeight_()
354 {
355  using Teuchos::as;
356 
357  int error = 0;
358  char title[100];
359  int exoid = 0;
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);
364 
365  if (typestr.compare("region") == 0) {
366  if (dimension_ == 3)
367  this->setEntityTypes(typestr, "vertex", "vertex");
368  else
369  // automatically downgrade primary entity if problem is only 2D
370  this->setEntityTypes("face", "vertex", "vertex");
371  }
372  else if (typestr.compare("vertex") == 0) {
373  if (dimension_ == 3)
374  this->setEntityTypes(typestr, "region", "region");
375  else
376  this->setEntityTypes(typestr, "face", "face");
377  }
378  else {
380  }
381 
382  coords_ = new scalar_t [num_nodes_ * dimension_];
383 
384  error += im_ex_get_coord(exoid, coords_, coords_ + num_nodes_,
385  coords_ + 2 * num_nodes_);
386 
387  element_num_map_ = new gno_t[num_elem_];
388  std::vector<int> tmp;
389  tmp.resize(num_elem_);
390 
391  // BDD cast to int did not always work!
392  // error += im_ex_get_elem_num_map(exoid, (int *)element_num_map_)
393  // This may be a case of calling the wrong method
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]);
397 
398  tmp.clear();
399  tmp.resize(num_nodes_);
400  node_num_map_ = new gno_t [num_nodes_];
401 
402  // BDD cast to int did not always work!
403  // error += im_ex_get_node_num_map(exoid, (int *)node_num_map_);
404  // This may be a case of calling the wrong method
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]);
408 
409  nodeTopology = new enum EntityTopologyType[num_nodes_];
410  for (int i=0;i<num_nodes_;i++)
411  nodeTopology[i] = POINT;
412  elemTopology = new enum EntityTopologyType[num_elem_];
413  for (int i=0;i<num_elem_;i++) {
414  if (dimension_==2)
415  elemTopology[i] = QUADRILATERAL;
416  else
417  elemTopology[i] = HEXAHEDRON;
418  }
419 
420  int *elem_blk_ids = new int [num_elem_blk];
421  error += im_ex_get_elem_blk_ids(exoid, elem_blk_ids);
422 
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];
428 
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];
436  }
437 
438  delete[] elem_type;
439  elem_type = NULL;
440  delete[] num_attr;
441  num_attr = NULL;
442  Acoords_ = new scalar_t [num_elem_ * dimension_];
443  int a = 0;
444  std::vector<std::vector<gno_t> > sur_elem;
445  sur_elem.resize(num_nodes_);
446 
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]);
450 
451  for(int i = 0; i < num_elem_this_blk[b]; i++) {
452  Acoords_[a] = 0;
453  Acoords_[num_elem_ + a] = 0;
454 
455  if (3 == dimension_) {
456  Acoords_[2 * num_elem_ + a] = 0;
457  }
458 
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];
463 
464  if(3 == dimension_) {
465  Acoords_[2 * num_elem_ + a] += coords_[2 * num_nodes_ + node];
466  }
467 
468  /*
469  * in the case of degenerate elements, where a node can be
470  * entered into the connect table twice, need to check to
471  * make sure that this element is not already listed as
472  * surrounding this node
473  */
474  if (sur_elem[node].empty() ||
475  element_num_map_[a] != sur_elem[node][sur_elem[node].size()-1]) {
476  /* Add the element to the list */
477  sur_elem[node].push_back(element_num_map_[a]);
478  }
479  }
480 
481  Acoords_[a] /= num_nodes_per_elem[b];
482  Acoords_[num_elem_ + a] /= num_nodes_per_elem[b];
483 
484  if(3 == dimension_) {
485  Acoords_[2 * num_elem_ + a] /= num_nodes_per_elem[b];
486  }
487 
488  a++;
489  }
490 
491  }
492 
493  delete[] elem_blk_ids;
494  elem_blk_ids = NULL;
495  int nnodes_per_elem = num_nodes_per_elem[0];
496  elemToNode_ = new gno_t[num_elem_ * nnodes_per_elem];
497  int telct = 0;
498  elemOffsets_ = new lno_t [num_elem_+1];
499  tnoct_ = 0;
500  int **reconnect = new int * [num_elem_];
501  size_t max_nsur = 0;
502 
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]];
507 
508  for (int j = 0; j < num_nodes_per_elem[b]; j++) {
509  elemToNode_[tnoct_]=
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];
512  ++tnoct_;
513  }
514 
515  ++telct;
516  }
517  }
518  elemOffsets_[telct] = tnoct_;
519 
520  delete[] num_nodes_per_elem;
521  num_nodes_per_elem = NULL;
522  delete[] num_elem_this_blk;
523  num_elem_this_blk = NULL;
524 
525  for(int b = 0; b < num_elem_blk; b++) {
526  delete[] connect[b];
527  }
528 
529  delete[] connect;
530  connect = NULL;
531 
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];
535 
536  /* Allocate memory necessary for the adjacency */
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;
541 
542  for (int i=0; i < max_side_nodes; i++) {
543  side_nodes[i]=-999;
544  mirror_nodes[i]=-999;
545  }
546 
547  /* Find the adjacency for a nodal based decomposition */
548  nEadj_ = 0;
549  nNadj_ = 0;
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"
553  << std::endl;
554  } else {
555  size_t nsur = sur_elem[ncnt].size();
556  if (nsur > max_nsur)
557  max_nsur = nsur;
558  }
559  }
560 
561  nodeToElem_ = new gno_t[num_nodes_ * max_nsur];
562  nodeOffsets_ = new lno_t[num_nodes_+1];
563  telct_ = 0;
564 
565  for (int ncnt = 0; ncnt < num_nodes_; ncnt++) {
566  nodeOffsets_[ncnt] = telct_;
567  nStart_[ncnt] = nNadj_;
568  MapType nAdjMap;
569 
570  for (size_t i = 0; i < sur_elem[ncnt].size(); i++) {
571  nodeToElem_[telct_] = sur_elem[ncnt][i];
572  ++telct_;
573 
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]);
580 
581  if (node_num_map_[ncnt] != elemToNode_[elemOffsets_[ecnt]+j] &&
582  iter == nAdjMap.end() ) {
583  nAdj.push_back(elemToNode_[elemOffsets_[ecnt]+j]);
584  nNadj_++;
585  nAdjMap.insert({elemToNode_[elemOffsets_[ecnt]+j],
586  elemToNode_[elemOffsets_[ecnt]+j]});
587  }
588  }
589 
590  break;
591  }
592  }
593 #endif
594  }
595 
596  nAdjMap.clear();
597  }
598 
599  nodeOffsets_[num_nodes_] = telct_;
600  nStart_[num_nodes_] = nNadj_;
601 
602  nAdj_ = new gno_t [nNadj_];
603 
604  for (size_t i=0; i < nNadj_; i++) {
605  nAdj_[i] = as<gno_t>(nAdj[i]);
606  }
607 
608  int nprocs = comm.getSize();
609  //if (nprocs > 1) {
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);
614 
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;
617  int proc = 0;
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);
622 
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;
633 
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);
641  }
642  delete[] node_cmap_ids;
643  node_cmap_ids = NULL;
644  int *sendCount = new int [nprocs];
645  int *recvCount = new int [nprocs];
646 
647  // Post receives
648  RCP<CommRequest<int> >*requests=new RCP<CommRequest<int> >[num_node_cmaps];
649  for (int cnt = 0, i = 0; i < num_node_cmaps; i++) {
650  try {
651  requests[cnt++] =
652  Teuchos::ireceive<int,int>(comm,
653  rcp(&(recvCount[node_proc_ids[i][0]]),
654  false),
655  node_proc_ids[i][0]);
656  }
658  }
659 
660  Teuchos::barrier<int>(comm);
661  size_t totalsend = 0;
662 
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;
667  }
668  totalsend += sendCount[node_proc_ids[j][0]];
669  }
670 
671  // Send data; can use readySend since receives are posted.
672  for (int i = 0; i < num_node_cmaps; i++) {
673  try {
674  Teuchos::readySend<int,int>(comm, sendCount[node_proc_ids[i][0]],
675  node_proc_ids[i][0]);
676  }
678  }
679 
680  // Wait for messages to return.
681  try {
682  Teuchos::waitAll<int>(comm, arrayView(requests, num_node_cmaps));
683  }
685 
686  delete [] requests;
687 
688  // Allocate the receive buffer.
689  size_t totalrecv = 0;
690  int maxMsg = 0;
691  int nrecvranks = 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]];
695  nrecvranks++;
696  if (recvCount[node_proc_ids[i][0]] > maxMsg)
697  maxMsg = recvCount[node_proc_ids[i][0]];
698  }
699  }
700 
701  gno_t *rbuf = NULL;
702  if (totalrecv) rbuf = new gno_t[totalrecv];
703 
704  requests = new RCP<CommRequest<int> > [nrecvranks];
705 
706  // Error checking for memory and message size.
707  int OK[2] = {1,1};
708  // OK[0] -- true/false indicating whether each message size fits in an int
709  // (for MPI).
710  // OK[1] -- true/false indicating whether memory allocs are OK
711  int gOK[2]; // For global reduce of OK.
712 
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;
716 
717  // Post receives
718 
719  size_t offset = 0;
720 
721  if (OK[0] && OK[1]) {
722  int rcnt = 0;
723  for (int i = 0; i < num_node_cmaps; i++) {
724  if (recvCount[node_proc_ids[i][0]]) {
725  try {
726  requests[rcnt++] =
727  Teuchos::
728  ireceive<int,gno_t>(comm,
729  Teuchos::arcp(&rbuf[offset], 0,
730  recvCount[node_proc_ids[i][0]],
731  false),
732  node_proc_ids[i][0]);
733  }
735  }
736  offset += recvCount[node_proc_ids[i][0]];
737  }
738  }
739 
740  delete[] recvCount;
741 
742  // Use barrier for error checking
743  Teuchos::reduceAll<int>(comm, Teuchos::REDUCE_MIN, 2, OK, gOK);
744  if (!gOK[0] || !gOK[1]) {
745  delete [] rbuf;
746  delete [] requests;
747  if (!gOK[0])
748  throw std::runtime_error("Max single message length exceeded");
749  else
750  throw std::bad_alloc();
751  }
752 
753  gno_t *sbuf = NULL;
754  if (totalsend) sbuf = new gno_t[totalsend];
755  a = 0;
756 
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];
764  }
765  }
766  }
767 
768  delete[] node_cmap_node_cnts;
769  node_cmap_node_cnts = NULL;
770 
771  for(int j = 0; j < num_node_cmaps; j++) {
772  delete[] node_ids[j];
773  }
774 
775  delete[] node_ids;
776  node_ids = NULL;
777  ArrayRCP<gno_t> sendBuf;
778 
779  if (totalsend)
780  sendBuf = ArrayRCP<gno_t>(sbuf, 0, totalsend, true);
781  else
782  sendBuf = Teuchos::null;
783 
784  // Send data; can use readySend since receives are posted.
785  offset = 0;
786  for (int i = 0; i < num_node_cmaps; i++) {
787  if (sendCount[node_proc_ids[i][0]]) {
788  try{
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]);
793  }
795  }
796  offset += sendCount[node_proc_ids[i][0]];
797  }
798 
799  for(int j = 0; j < num_node_cmaps; j++) {
800  delete[] node_proc_ids[j];
801  }
802 
803  delete[] node_proc_ids;
804  node_proc_ids = NULL;
805  delete[] sendCount;
806 
807  // Wait for messages to return.
808  try{
809  Teuchos::waitAll<int>(comm, Teuchos::arrayView(requests, nrecvranks));
810  }
812 
813  delete[] requests;
814  a = 0;
815 
816  for (int i = 0; i < num_node_cmaps; i++) {
817  int num_nodes_this_processor = rbuf[a++];
818 
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++];
822 
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++]);
827  }
828 
829  break;
830  }
831  }
832  }
833  }
834 
835  delete[] rbuf;
836  //}
837 
838 #ifndef USE_MESH_ADAPTER
839  for(int ecnt=0; ecnt < num_elem_; ecnt++) {
840  eStart_[ecnt] = nEadj_;
841  MapType eAdjMap;
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);
848 
849  if(element_num_map_[ecnt] != entry &&
850  iter == eAdjMap.end()) {
851  eAdj.push_back(entry);
852  nEadj_++;
853  eAdjMap.insert({entry, entry});
854  }
855  }
856  }
857 
858  eAdjMap.clear();
859  }
860 #endif
861 
862  for(int b = 0; b < num_elem_; b++) {
863  delete[] reconnect[b];
864  }
865 
866  delete[] reconnect;
867  reconnect = NULL;
868  eStart_[num_elem_] = nEadj_;
869 
870  eAdj_ = new gno_t [nEadj_];
871 
872  for (size_t i=0; i < nEadj_; i++) {
873  eAdj_[i] = as<gno_t>(eAdj[i]);
874  }
875 
876  delete[] side_nodes;
877  side_nodes = NULL;
878  delete[] mirror_nodes;
879  mirror_nodes = NULL;
880 
881  if (nWeightsPerEntity_ > 0) {
882  entityDegreeWeight_ = new bool [nWeightsPerEntity_];
883  for (int i=0; i < nWeightsPerEntity_; i++) {
884  entityDegreeWeight_[i] = false;
885  }
886  }
887 }
888 
890 template <typename User>
892 {
893  if (idx >= 0 && idx < nWeightsPerEntity_)
894  entityDegreeWeight_[idx] = true;
895  else
896  std::cout << "WARNING: invalid entity weight index, " << idx << ", ignored"
897  << std::endl;
898 }
899 
900 template <typename User>
902 {
903  std::string fn(" PamgenMesh ");
904  std::cout << me << fn
905  << " dim = " << dimension_
906  << " nodes = " << num_nodes_
907  << " nelems = " << num_elem_
908  << std::endl;
909 
910  for (int i = 0; i < num_elem_; i++) {
911  std::cout << me << fn << i
912  << " Elem " << element_num_map_[i]
913  << " Coords: ";
914  for (int j = 0; j < dimension_; j++)
915  std::cout << Acoords_[i + j * num_elem_] << " ";
916  std::cout << std::endl;
917  }
918 
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]
923  << " Graph: ";
924  for (int j = eStart_[i]; j < eStart_[i+1]; j++)
925  std::cout << eAdj_[j] << " ";
926  std::cout << std::endl;
927  }
928 #endif
929 }
930 
931 } //namespace Zoltan2
932 
933 #endif
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.
default_part_t part_t
The data type to represent part numbers.
void getTopologyViewOf(MeshEntityType etype, enum EntityTopologyType const *&Types) const
Provide a pointer to the entity topology types.
#define Z2_THROW_NOT_IMPLEMENTED
static ArrayRCP< ArrayRCP< zscalar_t > > weights
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 ...
default_lno_t lno_t
The ordinal type (e.g., int, long, int64_t) that represents local counts and local indices...
void get2ndAdjsView(MeshEntityType sourcetarget, MeshEntityType through, const lno_t *&offsets, const gno_t *&adjacencyIds) const
if avail2ndAdjs(), set pointers to this process&#39; second adjacencies
void getIDsViewOf(MeshEntityType etype, const gno_t *&Ids) const
Provide a pointer to this process&#39; identifiers.
EntityTopologyType
Enumerate entity topology types for meshes: points,lines,polygons,triangles,quadrilaterals, polyhedrons, tetrahedrons, hexhedrons, prisms, or pyramids.
default_gno_t gno_t
The ordinal type (e.g., int, long, int64_t) that can represent global counts and identifiers.
InputTraits< User >::node_t node_t
default_node_t node_t
The Kokkos node type. This is only meaningful for users of Tpetra objects.
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&#39; 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&#39; mesh first adjacencies.
This class represents a mesh.
default_scalar_t scalar_t
The data type for weights and coordinates.
This file defines the StridedData class.