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;
101  typedef User user_t;
102  typedef std::map<gno_t, gno_t> MapType;
103 
111  PamgenMeshAdapter(const Comm<int> &comm, std::string typestr="region",
112  int nEntWgts=0);
113 
120  void setWeightIsDegree(int idx);
121 
122  void print(int);
123 
125  // The MeshAdapter interface.
126  // This is the interface that would be called by a model or a problem .
128 
129  bool areEntityIDsUnique(MeshEntityType etype) const {
130  return etype==MESH_REGION;
131  }
132 
133  size_t getLocalNumOf(MeshEntityType etype) const
134  {
135  if ((MESH_REGION == etype && 3 == dimension_) ||
136  (MESH_FACE == etype && 2 == dimension_)) {
137  return num_elem_;
138  }
139 
140  if (MESH_VERTEX == etype) {
141  return num_nodes_;
142  }
143 
144  return 0;
145  }
146 
147  void getIDsViewOf(MeshEntityType etype, const gno_t *&Ids) const
148  {
149  if ((MESH_REGION == etype && 3 == dimension_) ||
150  (MESH_FACE == etype && 2 == dimension_)) {
151  Ids = element_num_map_;
152  }
153 
154  else if (MESH_VERTEX == etype) {
155  Ids = node_num_map_;
156  }
157 
158  else Ids = NULL;
159  }
160 
162  enum EntityTopologyType const *&Types) const {
163  if ((MESH_REGION == etype && 3 == dimension_) ||
164  (MESH_FACE == etype && 2 == dimension_)) {
165  Types = elemTopology;
166  }
167 
168  else if (MESH_VERTEX == etype) {
169  Types = nodeTopology;
170  }
171 
172  else Types = NULL;
173  }
174 
175  void getWeightsViewOf(MeshEntityType etype, const scalar_t *&weights,
176  int &stride, int idx = 0) const
177  {
178  weights = NULL;
179  stride = 0;
180  }
181 
182  int getDimension() const { return dimension_; }
183 
184  void getCoordinatesViewOf(MeshEntityType etype, const scalar_t *&coords,
185  int &stride, int dim) const {
186  if ((MESH_REGION == etype && 3 == dimension_) ||
187  (MESH_FACE == etype && 2 == dimension_)) {
188  if (dim == 0) {
189  coords = Acoords_;
190  } else if (dim == 1) {
191  coords = Acoords_ + num_elem_;
192  } else if (dim == 2) {
193  coords = Acoords_ + 2 * num_elem_;
194  }
195  stride = 1;
196  } else if (MESH_REGION == etype && 2 == dimension_) {
197  coords = NULL;
198  stride = 0;
199  } else if (MESH_VERTEX == etype) {
200  if (dim == 0) {
201  coords = coords_;
202  } else if (dim == 1) {
203  coords = coords_ + num_nodes_;
204  } else if (dim == 2) {
205  coords = coords_ + 2 * num_nodes_;
206  }
207  stride = 1;
208  } else {
209  coords = NULL;
210  stride = 0;
212  }
213  }
214 
215  bool availAdjs(MeshEntityType source, MeshEntityType target) const {
216  if ((MESH_REGION == source && MESH_VERTEX == target && 3 == dimension_) ||
217  (MESH_FACE == source && MESH_VERTEX == target && 2 == dimension_) ||
218  (MESH_VERTEX == source && MESH_REGION == target && 3 == dimension_) ||
219  (MESH_VERTEX == source && MESH_FACE == target && 2 == dimension_)) {
220  return TRUE;
221  }
222 
223  return false;
224  }
225 
226  size_t getLocalNumAdjs(MeshEntityType source, MeshEntityType target) const
227  {
228  if ((MESH_REGION == source && MESH_VERTEX == target && 3 == dimension_) ||
229  (MESH_FACE == source && MESH_VERTEX == target && 2 == dimension_)) {
230  return tnoct_;
231  }
232 
233  if ((MESH_VERTEX == source && MESH_REGION == target && 3 == dimension_) ||
234  (MESH_VERTEX == source && MESH_FACE == target && 2 == dimension_)) {
235  return telct_;
236  }
237 
238  return 0;
239  }
240 
242  const lno_t *&offsets, const gno_t *& adjacencyIds) const
243  {
244  if ((MESH_REGION == source && MESH_VERTEX == target && 3 == dimension_) ||
245  (MESH_FACE == source && MESH_VERTEX == target && 2 == dimension_)) {
246  offsets = elemOffsets_;
247  adjacencyIds = elemToNode_;
248  } else if ((MESH_REGION==target && MESH_VERTEX==source && 3==dimension_) ||
249  (MESH_FACE==target && MESH_VERTEX==source && 2==dimension_)) {
250  offsets = nodeOffsets_;
251  adjacencyIds = nodeToElem_;
252  } else if (MESH_REGION == source && 2 == dimension_) {
253  offsets = NULL;
254  adjacencyIds = NULL;
255  } else {
256  offsets = NULL;
257  adjacencyIds = NULL;
259  }
260  }
261 
262  //#define USE_MESH_ADAPTER
263 #ifndef USE_MESH_ADAPTER
264  bool avail2ndAdjs(MeshEntityType sourcetarget, MeshEntityType through) const
265  {
266  if (through == MESH_VERTEX) {
267  if (sourcetarget == MESH_REGION && dimension_ == 3) return true;
268  if (sourcetarget == MESH_FACE && dimension_ == 2) return true;
269  }
270  if (sourcetarget == MESH_VERTEX) {
271  if (through == MESH_REGION && dimension_ == 3) return true;
272  if (through == MESH_FACE && dimension_ == 2) return true;
273  }
274  return false;
275  }
276 
277  size_t getLocalNum2ndAdjs(MeshEntityType sourcetarget,
278  MeshEntityType through) const
279  {
280  if (through == MESH_VERTEX &&
281  ((sourcetarget == MESH_REGION && dimension_ == 3) ||
282  (sourcetarget == MESH_FACE && dimension_ == 2))) {
283  return nEadj_;
284  }
285 
286  if (sourcetarget == MESH_VERTEX &&
287  ((through == MESH_REGION && dimension_ == 3) ||
288  (through == MESH_FACE && dimension_ == 2))) {
289  return nNadj_;
290  }
291 
292  return 0;
293  }
294 
295  void get2ndAdjsView(MeshEntityType sourcetarget, MeshEntityType through,
296  const lno_t *&offsets, const gno_t *&adjacencyIds) const
297  {
298  if (through == MESH_VERTEX &&
299  ((sourcetarget == MESH_REGION && dimension_ == 3) ||
300  (sourcetarget == MESH_FACE && dimension_ == 2))) {
301  offsets = eStart_;
302  adjacencyIds = eAdj_;
303  } else if (sourcetarget == MESH_VERTEX &&
304  ((through == MESH_REGION && dimension_ == 3) ||
305  (through == MESH_FACE && dimension_ == 2))) {
306  offsets = nStart_;
307  adjacencyIds = nAdj_;
308  } else {
309  offsets = NULL;
310  adjacencyIds = NULL;
312  }
313  }
314 #endif
315 
316  bool useDegreeAsWeightOf(MeshEntityType etype, int idx) const
317  {
318  if ((MESH_REGION == etype && 3 == dimension_) ||
319  (MESH_FACE == etype && 2 == dimension_) ||
320  (MESH_VERTEX == etype)) {
321  return entityDegreeWeight_[idx];
322  }
323 
324  return false;
325  }
326 
327 private:
328  int dimension_, num_nodes_global_, num_elems_global_, num_nodes_, num_elem_;
329  gno_t *element_num_map_, *node_num_map_;
330  gno_t *elemToNode_;
331  lno_t tnoct_, *elemOffsets_;
332  gno_t *nodeToElem_;
333  lno_t telct_, *nodeOffsets_;
334 
335  int nWeightsPerEntity_;
336  bool *entityDegreeWeight_;
337 
338  scalar_t *coords_, *Acoords_;
339  lno_t *eStart_, *nStart_;
340  gno_t *eAdj_, *nAdj_;
341  size_t nEadj_, nNadj_;
342  EntityTopologyType* nodeTopology;
343  EntityTopologyType* elemTopology;
344 
345 };
346 
348 // Definitions
350 
351 template <typename User>
353  std::string typestr, int nEntWgts):
354  dimension_(0), nWeightsPerEntity_(nEntWgts), entityDegreeWeight_()
355 {
356  using Teuchos::as;
357 
358  int error = 0;
359  char title[100];
360  int exoid = 0;
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);
365 
366  if (typestr.compare("region") == 0) {
367  if (dimension_ == 3)
368  this->setEntityTypes(typestr, "vertex", "vertex");
369  else
370  // automatically downgrade primary entity if problem is only 2D
371  this->setEntityTypes("face", "vertex", "vertex");
372  }
373  else if (typestr.compare("vertex") == 0) {
374  if (dimension_ == 3)
375  this->setEntityTypes(typestr, "region", "region");
376  else
377  this->setEntityTypes(typestr, "face", "face");
378  }
379  else {
381  }
382 
383  coords_ = new scalar_t [num_nodes_ * dimension_];
384 
385  error += im_ex_get_coord(exoid, coords_, coords_ + num_nodes_,
386  coords_ + 2 * num_nodes_);
387 
388  element_num_map_ = new gno_t[num_elem_];
389  std::vector<int> tmp;
390  tmp.resize(num_elem_);
391 
392  // BDD cast to int did not always work!
393  // error += im_ex_get_elem_num_map(exoid, (int *)element_num_map_)
394  // This may be a case of calling the wrong method
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]);
398 
399  tmp.clear();
400  tmp.resize(num_nodes_);
401  node_num_map_ = new gno_t [num_nodes_];
402 
403  // BDD cast to int did not always work!
404  // error += im_ex_get_node_num_map(exoid, (int *)node_num_map_);
405  // This may be a case of calling the wrong method
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]);
409 
410  nodeTopology = new enum EntityTopologyType[num_nodes_];
411  for (int i=0;i<num_nodes_;i++)
412  nodeTopology[i] = POINT;
413  elemTopology = new enum EntityTopologyType[num_elem_];
414  for (int i=0;i<num_elem_;i++) {
415  if (dimension_==2)
416  elemTopology[i] = QUADRILATERAL;
417  else
418  elemTopology[i] = HEXAHEDRON;
419  }
420 
421  int *elem_blk_ids = new int [num_elem_blk];
422  error += im_ex_get_elem_blk_ids(exoid, elem_blk_ids);
423 
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];
429 
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];
437  }
438 
439  delete[] elem_type;
440  elem_type = NULL;
441  delete[] num_attr;
442  num_attr = NULL;
443  Acoords_ = new scalar_t [num_elem_ * dimension_];
444  int a = 0;
445  std::vector<std::vector<gno_t> > sur_elem;
446  sur_elem.resize(num_nodes_);
447 
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]);
451 
452  for(int i = 0; i < num_elem_this_blk[b]; i++) {
453  Acoords_[a] = 0;
454  Acoords_[num_elem_ + a] = 0;
455 
456  if (3 == dimension_) {
457  Acoords_[2 * num_elem_ + a] = 0;
458  }
459 
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];
464 
465  if(3 == dimension_) {
466  Acoords_[2 * num_elem_ + a] += coords_[2 * num_nodes_ + node];
467  }
468 
469  /*
470  * in the case of degenerate elements, where a node can be
471  * entered into the connect table twice, need to check to
472  * make sure that this element is not already listed as
473  * surrounding this node
474  */
475  if (sur_elem[node].empty() ||
476  element_num_map_[a] != sur_elem[node][sur_elem[node].size()-1]) {
477  /* Add the element to the list */
478  sur_elem[node].push_back(element_num_map_[a]);
479  }
480  }
481 
482  Acoords_[a] /= num_nodes_per_elem[b];
483  Acoords_[num_elem_ + a] /= num_nodes_per_elem[b];
484 
485  if(3 == dimension_) {
486  Acoords_[2 * num_elem_ + a] /= num_nodes_per_elem[b];
487  }
488 
489  a++;
490  }
491 
492  }
493 
494  delete[] elem_blk_ids;
495  elem_blk_ids = NULL;
496  int nnodes_per_elem = num_nodes_per_elem[0];
497  elemToNode_ = new gno_t[num_elem_ * nnodes_per_elem];
498  int telct = 0;
499  elemOffsets_ = new lno_t [num_elem_+1];
500  tnoct_ = 0;
501  int **reconnect = new int * [num_elem_];
502  size_t max_nsur = 0;
503 
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]];
508 
509  for (int j = 0; j < num_nodes_per_elem[b]; j++) {
510  elemToNode_[tnoct_]=
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];
513  ++tnoct_;
514  }
515 
516  ++telct;
517  }
518  }
519  elemOffsets_[telct] = tnoct_;
520 
521  delete[] num_nodes_per_elem;
522  num_nodes_per_elem = NULL;
523  delete[] num_elem_this_blk;
524  num_elem_this_blk = NULL;
525 
526  for(int b = 0; b < num_elem_blk; b++) {
527  delete[] connect[b];
528  }
529 
530  delete[] connect;
531  connect = NULL;
532 
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];
536 
537  /* Allocate memory necessary for the adjacency */
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;
542 
543  for (int i=0; i < max_side_nodes; i++) {
544  side_nodes[i]=-999;
545  mirror_nodes[i]=-999;
546  }
547 
548  /* Find the adjacency for a nodal based decomposition */
549  nEadj_ = 0;
550  nNadj_ = 0;
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"
554  << std::endl;
555  } else {
556  size_t nsur = sur_elem[ncnt].size();
557  if (nsur > max_nsur)
558  max_nsur = nsur;
559  }
560  }
561 
562  nodeToElem_ = new gno_t[num_nodes_ * max_nsur];
563  nodeOffsets_ = new lno_t[num_nodes_+1];
564  telct_ = 0;
565 
566  for (int ncnt = 0; ncnt < num_nodes_; ncnt++) {
567  nodeOffsets_[ncnt] = telct_;
568  nStart_[ncnt] = nNadj_;
569  MapType nAdjMap;
570 
571  for (size_t i = 0; i < sur_elem[ncnt].size(); i++) {
572  nodeToElem_[telct_] = sur_elem[ncnt][i];
573  ++telct_;
574 
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]);
581 
582  if (node_num_map_[ncnt] != elemToNode_[elemOffsets_[ecnt]+j] &&
583  iter == nAdjMap.end() ) {
584  nAdj.push_back(elemToNode_[elemOffsets_[ecnt]+j]);
585  nNadj_++;
586  nAdjMap.insert({elemToNode_[elemOffsets_[ecnt]+j],
587  elemToNode_[elemOffsets_[ecnt]+j]});
588  }
589  }
590 
591  break;
592  }
593  }
594 #endif
595  }
596 
597  nAdjMap.clear();
598  }
599 
600  nodeOffsets_[num_nodes_] = telct_;
601  nStart_[num_nodes_] = nNadj_;
602 
603  nAdj_ = new gno_t [nNadj_];
604 
605  for (size_t i=0; i < nNadj_; i++) {
606  nAdj_[i] = as<gno_t>(nAdj[i]);
607  }
608 
609  int nprocs = comm.getSize();
610  //if (nprocs > 1) {
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);
615 
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;
618  int proc = 0;
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);
623 
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;
634 
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);
642  }
643  delete[] node_cmap_ids;
644  node_cmap_ids = NULL;
645  int *sendCount = new int [nprocs];
646  int *recvCount = new int [nprocs];
647 
648  // Post receives
649  RCP<CommRequest<int> >*requests=new RCP<CommRequest<int> >[num_node_cmaps];
650  for (int cnt = 0, i = 0; i < num_node_cmaps; i++) {
651  try {
652  requests[cnt++] =
653  Teuchos::ireceive<int,int>(comm,
654  rcp(&(recvCount[node_proc_ids[i][0]]),
655  false),
656  node_proc_ids[i][0]);
657  }
659  }
660 
661  Teuchos::barrier<int>(comm);
662  size_t totalsend = 0;
663 
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;
668  }
669  totalsend += sendCount[node_proc_ids[j][0]];
670  }
671 
672  // Send data; can use readySend since receives are posted.
673  for (int i = 0; i < num_node_cmaps; i++) {
674  try {
675  Teuchos::readySend<int,int>(comm, sendCount[node_proc_ids[i][0]],
676  node_proc_ids[i][0]);
677  }
679  }
680 
681  // Wait for messages to return.
682  try {
683  Teuchos::waitAll<int>(comm, arrayView(requests, num_node_cmaps));
684  }
686 
687  delete [] requests;
688 
689  // Allocate the receive buffer.
690  size_t totalrecv = 0;
691  int maxMsg = 0;
692  int nrecvranks = 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]];
696  nrecvranks++;
697  if (recvCount[node_proc_ids[i][0]] > maxMsg)
698  maxMsg = recvCount[node_proc_ids[i][0]];
699  }
700  }
701 
702  gno_t *rbuf = NULL;
703  if (totalrecv) rbuf = new gno_t[totalrecv];
704 
705  requests = new RCP<CommRequest<int> > [nrecvranks];
706 
707  // Error checking for memory and message size.
708  int OK[2] = {1,1};
709  // OK[0] -- true/false indicating whether each message size fits in an int
710  // (for MPI).
711  // OK[1] -- true/false indicating whether memory allocs are OK
712  int gOK[2]; // For global reduce of OK.
713 
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;
717 
718  // Post receives
719 
720  size_t offset = 0;
721 
722  if (OK[0] && OK[1]) {
723  int rcnt = 0;
724  for (int i = 0; i < num_node_cmaps; i++) {
725  if (recvCount[node_proc_ids[i][0]]) {
726  try {
727  requests[rcnt++] =
728  Teuchos::
729  ireceive<int,gno_t>(comm,
730  Teuchos::arcp(&rbuf[offset], 0,
731  recvCount[node_proc_ids[i][0]],
732  false),
733  node_proc_ids[i][0]);
734  }
736  }
737  offset += recvCount[node_proc_ids[i][0]];
738  }
739  }
740 
741  delete[] recvCount;
742 
743  // Use barrier for error checking
744  Teuchos::reduceAll<int>(comm, Teuchos::REDUCE_MIN, 2, OK, gOK);
745  if (!gOK[0] || !gOK[1]) {
746  delete [] rbuf;
747  delete [] requests;
748  if (!gOK[0])
749  throw std::runtime_error("Max single message length exceeded");
750  else
751  throw std::bad_alloc();
752  }
753 
754  gno_t *sbuf = NULL;
755  if (totalsend) sbuf = new gno_t[totalsend];
756  a = 0;
757 
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];
765  }
766  }
767  }
768 
769  delete[] node_cmap_node_cnts;
770  node_cmap_node_cnts = NULL;
771 
772  for(int j = 0; j < num_node_cmaps; j++) {
773  delete[] node_ids[j];
774  }
775 
776  delete[] node_ids;
777  node_ids = NULL;
778  ArrayRCP<gno_t> sendBuf;
779 
780  if (totalsend)
781  sendBuf = ArrayRCP<gno_t>(sbuf, 0, totalsend, true);
782  else
783  sendBuf = Teuchos::null;
784 
785  // Send data; can use readySend since receives are posted.
786  offset = 0;
787  for (int i = 0; i < num_node_cmaps; i++) {
788  if (sendCount[node_proc_ids[i][0]]) {
789  try{
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]);
794  }
796  }
797  offset += sendCount[node_proc_ids[i][0]];
798  }
799 
800  for(int j = 0; j < num_node_cmaps; j++) {
801  delete[] node_proc_ids[j];
802  }
803 
804  delete[] node_proc_ids;
805  node_proc_ids = NULL;
806  delete[] sendCount;
807 
808  // Wait for messages to return.
809  try{
810  Teuchos::waitAll<int>(comm, Teuchos::arrayView(requests, nrecvranks));
811  }
813 
814  delete[] requests;
815  a = 0;
816 
817  for (int i = 0; i < num_node_cmaps; i++) {
818  int num_nodes_this_processor = rbuf[a++];
819 
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++];
823 
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++]);
828  }
829 
830  break;
831  }
832  }
833  }
834  }
835 
836  delete[] rbuf;
837  //}
838 
839 #ifndef USE_MESH_ADAPTER
840  for(int ecnt=0; ecnt < num_elem_; ecnt++) {
841  eStart_[ecnt] = nEadj_;
842  MapType eAdjMap;
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);
849 
850  if(element_num_map_[ecnt] != entry &&
851  iter == eAdjMap.end()) {
852  eAdj.push_back(entry);
853  nEadj_++;
854  eAdjMap.insert({entry, entry});
855  }
856  }
857  }
858 
859  eAdjMap.clear();
860  }
861 #endif
862 
863  for(int b = 0; b < num_elem_; b++) {
864  delete[] reconnect[b];
865  }
866 
867  delete[] reconnect;
868  reconnect = NULL;
869  eStart_[num_elem_] = nEadj_;
870 
871  eAdj_ = new gno_t [nEadj_];
872 
873  for (size_t i=0; i < nEadj_; i++) {
874  eAdj_[i] = as<gno_t>(eAdj[i]);
875  }
876 
877  delete[] side_nodes;
878  side_nodes = NULL;
879  delete[] mirror_nodes;
880  mirror_nodes = NULL;
881 
882  if (nWeightsPerEntity_ > 0) {
883  entityDegreeWeight_ = new bool [nWeightsPerEntity_];
884  for (int i=0; i < nWeightsPerEntity_; i++) {
885  entityDegreeWeight_[i] = false;
886  }
887  }
888 }
889 
891 template <typename User>
893 {
894  if (idx >= 0 && idx < nWeightsPerEntity_)
895  entityDegreeWeight_[idx] = true;
896  else
897  std::cout << "WARNING: invalid entity weight index, " << idx << ", ignored"
898  << std::endl;
899 }
900 
901 template <typename User>
903 {
904  std::string fn(" PamgenMesh ");
905  std::cout << me << fn
906  << " dim = " << dimension_
907  << " nodes = " << num_nodes_
908  << " nelems = " << num_elem_
909  << std::endl;
910 
911  for (int i = 0; i < num_elem_; i++) {
912  std::cout << me << fn << i
913  << " Elem " << element_num_map_[i]
914  << " Coords: ";
915  for (int j = 0; j < dimension_; j++)
916  std::cout << Acoords_[i + j * num_elem_] << " ";
917  std::cout << std::endl;
918  }
919 
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]
924  << " Graph: ";
925  for (int j = eStart_[i]; j < eStart_[i+1]; j++)
926  std::cout << eAdj_[j] << " ";
927  std::cout << std::endl;
928  }
929 #endif
930 }
931 
932 } //namespace Zoltan2
933 
934 #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...
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.
default_part_t part_t
The data type to represent part numbers.
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.
void getIDsViewOf(MeshEntityType etype, const gno_t *&Ids) const
Provide a pointer to this process&#39; identifiers.
static ArrayRCP< ArrayRCP< zscalar_t > > weights
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
default_lno_t lno_t
The ordinal type (e.g., int, long, int64_t) that represents local counts and local indices...
size_t getLocalNum2ndAdjs(MeshEntityType sourcetarget, MeshEntityType through) const
if avail2ndAdjs(), returns the number of second adjacencies on this process.
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 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.
default_scalar_t scalar_t
The data type for weights and coordinates.
bool availAdjs(MeshEntityType source, MeshEntityType target) const
Returns whether a first adjacency combination is available.
This file defines the StridedData class.