Zoltan2
Loading...
Searching...
No Matches
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
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
62namespace Zoltan2 {
63
90template <typename User>
91 class PamgenMeshAdapter: public MeshAdapter<User> {
92
93public:
94
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
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
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 offset_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 offset_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
327private:
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 offset_t tnoct_, *elemOffsets_;
332 gno_t *nodeToElem_;
333 offset_t telct_, *nodeOffsets_;
334
335 int nWeightsPerEntity_;
336 bool *entityDegreeWeight_;
337
338 scalar_t *coords_, *Acoords_;
339 offset_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
351template <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 double *dcoords = new double [num_nodes_ * dimension_];
384 coords_ = new scalar_t [num_nodes_ * dimension_];
385
386 error += im_ex_get_coord(exoid, dcoords, dcoords + num_nodes_,
387 dcoords + 2 * num_nodes_);
388
389 size_t dlen = num_nodes_ * dimension_;
390 for (size_t i = 0; i < dlen; i++) coords_[i] = as<scalar_t>(dcoords[i]);
391 delete [] dcoords;
392
393 element_num_map_ = new gno_t[num_elem_];
394 std::vector<int> tmp;
395 tmp.resize(num_elem_);
396
397 // BDD cast to int did not always work!
398 // error += im_ex_get_elem_num_map(exoid, (int *)element_num_map_)
399 // This may be a case of calling the wrong method
400 error += im_ex_get_elem_num_map(exoid, &tmp[0]);
401 for(size_t i = 0; i < tmp.size(); i++)
402 element_num_map_[i] = static_cast<gno_t>(tmp[i]);
403
404 tmp.clear();
405 tmp.resize(num_nodes_);
406 node_num_map_ = new gno_t [num_nodes_];
407
408 // BDD cast to int did not always work!
409 // error += im_ex_get_node_num_map(exoid, (int *)node_num_map_);
410 // This may be a case of calling the wrong method
411 error += im_ex_get_node_num_map(exoid, &tmp[0]);
412 for(size_t i = 0; i < tmp.size(); i++)
413 node_num_map_[i] = static_cast<gno_t>(tmp[i]);
414
415 nodeTopology = new enum EntityTopologyType[num_nodes_];
416 for (int i=0;i<num_nodes_;i++)
417 nodeTopology[i] = POINT;
418 elemTopology = new enum EntityTopologyType[num_elem_];
419 for (int i=0;i<num_elem_;i++) {
420 if (dimension_==2)
421 elemTopology[i] = QUADRILATERAL;
422 else
423 elemTopology[i] = HEXAHEDRON;
424 }
425
426 int *elem_blk_ids = new int [num_elem_blk];
427 error += im_ex_get_elem_blk_ids(exoid, elem_blk_ids);
428
429 int *num_nodes_per_elem = new int [num_elem_blk];
430 int *num_attr = new int [num_elem_blk];
431 int *num_elem_this_blk = new int [num_elem_blk];
432 char **elem_type = new char * [num_elem_blk];
433 int **connect = new int * [num_elem_blk];
434
435 for(int i = 0; i < num_elem_blk; i++){
436 elem_type[i] = new char [MAX_STR_LENGTH + 1];
437 error += im_ex_get_elem_block(exoid, elem_blk_ids[i], elem_type[i],
438 (int*)&(num_elem_this_blk[i]),
439 (int*)&(num_nodes_per_elem[i]),
440 (int*)&(num_attr[i]));
441 delete[] elem_type[i];
442 }
443
444 delete[] elem_type;
445 elem_type = NULL;
446 delete[] num_attr;
447 num_attr = NULL;
448 Acoords_ = new scalar_t [num_elem_ * dimension_];
449 int a = 0;
450 std::vector<std::vector<gno_t> > sur_elem;
451 sur_elem.resize(num_nodes_);
452
453 for(int b = 0; b < num_elem_blk; b++) {
454 connect[b] = new int [num_nodes_per_elem[b]*num_elem_this_blk[b]];
455 error += im_ex_get_elem_conn(exoid, elem_blk_ids[b], connect[b]);
456
457 for(int i = 0; i < num_elem_this_blk[b]; i++) {
458 Acoords_[a] = 0;
459 Acoords_[num_elem_ + a] = 0;
460
461 if (3 == dimension_) {
462 Acoords_[2 * num_elem_ + a] = 0;
463 }
464
465 for(int j = 0; j < num_nodes_per_elem[b]; j++) {
466 int node = connect[b][i * num_nodes_per_elem[b] + j] - 1;
467 Acoords_[a] += coords_[node];
468 Acoords_[num_elem_ + a] += coords_[num_nodes_ + node];
469
470 if(3 == dimension_) {
471 Acoords_[2 * num_elem_ + a] += coords_[2 * num_nodes_ + node];
472 }
473
474 /*
475 * in the case of degenerate elements, where a node can be
476 * entered into the connect table twice, need to check to
477 * make sure that this element is not already listed as
478 * surrounding this node
479 */
480 if (sur_elem[node].empty() ||
481 element_num_map_[a] != sur_elem[node][sur_elem[node].size()-1]) {
482 /* Add the element to the list */
483 sur_elem[node].push_back(element_num_map_[a]);
484 }
485 }
486
487 Acoords_[a] /= num_nodes_per_elem[b];
488 Acoords_[num_elem_ + a] /= num_nodes_per_elem[b];
489
490 if(3 == dimension_) {
491 Acoords_[2 * num_elem_ + a] /= num_nodes_per_elem[b];
492 }
493
494 a++;
495 }
496
497 }
498
499 delete[] elem_blk_ids;
500 elem_blk_ids = NULL;
501 int nnodes_per_elem = num_nodes_per_elem[0];
502 elemToNode_ = new gno_t[num_elem_ * nnodes_per_elem];
503 int telct = 0;
504 elemOffsets_ = new offset_t [num_elem_+1];
505 tnoct_ = 0;
506 int **reconnect = new int * [num_elem_];
507 size_t max_nsur = 0;
508
509 for (int b = 0; b < num_elem_blk; b++) {
510 for (int i = 0; i < num_elem_this_blk[b]; i++) {
511 elemOffsets_[telct] = tnoct_;
512 reconnect[telct] = new int [num_nodes_per_elem[b]];
513
514 for (int j = 0; j < num_nodes_per_elem[b]; j++) {
515 elemToNode_[tnoct_]=
516 as<gno_t>(node_num_map_[connect[b][i*num_nodes_per_elem[b] + j]-1]);
517 reconnect[telct][j] = connect[b][i*num_nodes_per_elem[b] + j];
518 ++tnoct_;
519 }
520
521 ++telct;
522 }
523 }
524 elemOffsets_[telct] = tnoct_;
525
526 delete[] num_nodes_per_elem;
527 num_nodes_per_elem = NULL;
528 delete[] num_elem_this_blk;
529 num_elem_this_blk = NULL;
530
531 for(int b = 0; b < num_elem_blk; b++) {
532 delete[] connect[b];
533 }
534
535 delete[] connect;
536 connect = NULL;
537
538 int max_side_nodes = nnodes_per_elem;
539 int *side_nodes = new int [max_side_nodes];
540 int *mirror_nodes = new int [max_side_nodes];
541
542 /* Allocate memory necessary for the adjacency */
543 eStart_ = new lno_t [num_elem_+1];
544 nStart_ = new lno_t [num_nodes_+1];
545 std::vector<int> eAdj;
546 std::vector<int> nAdj;
547
548 for (int i=0; i < max_side_nodes; i++) {
549 side_nodes[i]=-999;
550 mirror_nodes[i]=-999;
551 }
552
553 /* Find the adjacency for a nodal based decomposition */
554 nEadj_ = 0;
555 nNadj_ = 0;
556 for(int ncnt=0; ncnt < num_nodes_; ncnt++) {
557 if(sur_elem[ncnt].empty()) {
558 std::cout << "WARNING: Node = " << ncnt+1 << " has no elements"
559 << std::endl;
560 } else {
561 size_t nsur = sur_elem[ncnt].size();
562 if (nsur > max_nsur)
563 max_nsur = nsur;
564 }
565 }
566
567 nodeToElem_ = new gno_t[num_nodes_ * max_nsur];
568 nodeOffsets_ = new offset_t[num_nodes_+1];
569 telct_ = 0;
570
571 for (int ncnt = 0; ncnt < num_nodes_; ncnt++) {
572 nodeOffsets_[ncnt] = telct_;
573 nStart_[ncnt] = nNadj_;
574 MapType nAdjMap;
575
576 for (size_t i = 0; i < sur_elem[ncnt].size(); i++) {
577 nodeToElem_[telct_] = sur_elem[ncnt][i];
578 ++telct_;
579
580#ifndef USE_MESH_ADAPTER
581 for(int ecnt = 0; ecnt < num_elem_; ecnt++) {
582 if (element_num_map_[ecnt] == sur_elem[ncnt][i]) {
583 for (int j = 0; j < nnodes_per_elem; j++) {
584 typename MapType::iterator iter =
585 nAdjMap.find(elemToNode_[elemOffsets_[ecnt]+j]);
586
587 if (node_num_map_[ncnt] != elemToNode_[elemOffsets_[ecnt]+j] &&
588 iter == nAdjMap.end() ) {
589 nAdj.push_back(elemToNode_[elemOffsets_[ecnt]+j]);
590 nNadj_++;
591 nAdjMap.insert({elemToNode_[elemOffsets_[ecnt]+j],
592 elemToNode_[elemOffsets_[ecnt]+j]});
593 }
594 }
595
596 break;
597 }
598 }
599#endif
600 }
601
602 nAdjMap.clear();
603 }
604
605 nodeOffsets_[num_nodes_] = telct_;
606 nStart_[num_nodes_] = nNadj_;
607
608 nAdj_ = new gno_t [nNadj_];
609
610 for (size_t i=0; i < nNadj_; i++) {
611 nAdj_[i] = as<gno_t>(nAdj[i]);
612 }
613
614 int nprocs = comm.getSize();
615 //if (nprocs > 1) {
616 int neid=0,num_elem_blks_global,num_node_sets_global,num_side_sets_global;
617 error += im_ne_get_init_global(neid,&num_nodes_global_,&num_elems_global_,
618 &num_elem_blks_global,&num_node_sets_global,
619 &num_side_sets_global);
620
621 int num_internal_nodes, num_border_nodes, num_external_nodes;
622 int num_internal_elems, num_border_elems, num_node_cmaps, num_elem_cmaps;
623 int proc = 0;
624 error += im_ne_get_loadbal_param(neid, &num_internal_nodes,
625 &num_border_nodes, &num_external_nodes,
626 &num_internal_elems, &num_border_elems,
627 &num_node_cmaps, &num_elem_cmaps, proc);
628
629 int *node_cmap_ids = new int [num_node_cmaps];
630 int *node_cmap_node_cnts = new int [num_node_cmaps];
631 int *elem_cmap_ids = new int [num_elem_cmaps];
632 int *elem_cmap_elem_cnts = new int [num_elem_cmaps];
633 error += im_ne_get_cmap_params(neid, node_cmap_ids, node_cmap_node_cnts,
634 elem_cmap_ids, elem_cmap_elem_cnts, proc);
635 delete[] elem_cmap_ids;
636 elem_cmap_ids = NULL;
637 delete[] elem_cmap_elem_cnts;
638 elem_cmap_elem_cnts = NULL;
639
640 int **node_ids = new int * [num_node_cmaps];
641 int **node_proc_ids = new int * [num_node_cmaps];
642 for(int j = 0; j < num_node_cmaps; j++) {
643 node_ids[j] = new int [node_cmap_node_cnts[j]];
644 node_proc_ids[j] = new int [node_cmap_node_cnts[j]];
645 error += im_ne_get_node_cmap(neid, node_cmap_ids[j], node_ids[j],
646 node_proc_ids[j], proc);
647 }
648 delete[] node_cmap_ids;
649 node_cmap_ids = NULL;
650 int *sendCount = new int [nprocs];
651 int *recvCount = new int [nprocs];
652
653 // Post receives
654 RCP<CommRequest<int> >*requests=new RCP<CommRequest<int> >[num_node_cmaps];
655 for (int cnt = 0, i = 0; i < num_node_cmaps; i++) {
656 try {
657 requests[cnt++] =
658 Teuchos::ireceive<int,int>(comm,
659 rcp(&(recvCount[node_proc_ids[i][0]]),
660 false),
661 node_proc_ids[i][0]);
662 }
664 }
665
666 Teuchos::barrier<int>(comm);
667 size_t totalsend = 0;
668
669 for(int j = 0; j < num_node_cmaps; j++) {
670 sendCount[node_proc_ids[j][0]] = 1;
671 for(int i = 0; i < node_cmap_node_cnts[j]; i++) {
672 sendCount[node_proc_ids[j][i]] += sur_elem[node_ids[j][i]-1].size()+2;
673 }
674 totalsend += sendCount[node_proc_ids[j][0]];
675 }
676
677 // Send data; can use readySend since receives are posted.
678 for (int i = 0; i < num_node_cmaps; i++) {
679 try {
680 Teuchos::readySend<int,int>(comm, sendCount[node_proc_ids[i][0]],
681 node_proc_ids[i][0]);
682 }
684 }
685
686 // Wait for messages to return.
687 try {
688 Teuchos::waitAll<int>(comm, arrayView(requests, num_node_cmaps));
689 }
691
692 delete [] requests;
693
694 // Allocate the receive buffer.
695 size_t totalrecv = 0;
696 int maxMsg = 0;
697 int nrecvranks = 0;
698 for(int i = 0; i < num_node_cmaps; i++) {
699 if (recvCount[node_proc_ids[i][0]] > 0) {
700 totalrecv += recvCount[node_proc_ids[i][0]];
701 nrecvranks++;
702 if (recvCount[node_proc_ids[i][0]] > maxMsg)
703 maxMsg = recvCount[node_proc_ids[i][0]];
704 }
705 }
706
707 gno_t *rbuf = NULL;
708 if (totalrecv) rbuf = new gno_t[totalrecv];
709
710 requests = new RCP<CommRequest<int> > [nrecvranks];
711
712 // Error checking for memory and message size.
713 int OK[2] = {1,1};
714 // OK[0] -- true/false indicating whether each message size fits in an int
715 // (for MPI).
716 // OK[1] -- true/false indicating whether memory allocs are OK
717 int gOK[2]; // For global reduce of OK.
718
719 if (size_t(maxMsg) * sizeof(int) > INT_MAX) OK[0] = false;
720 if (totalrecv && !rbuf) OK[1] = 0;
721 if (!requests) OK[1] = 0;
722
723 // Post receives
724
725 size_t offset = 0;
726
727 if (OK[0] && OK[1]) {
728 int rcnt = 0;
729 for (int i = 0; i < num_node_cmaps; i++) {
730 if (recvCount[node_proc_ids[i][0]]) {
731 try {
732 requests[rcnt++] =
733 Teuchos::
734 ireceive<int,gno_t>(comm,
735 Teuchos::arcp(&rbuf[offset], 0,
736 recvCount[node_proc_ids[i][0]],
737 false),
738 node_proc_ids[i][0]);
739 }
741 }
742 offset += recvCount[node_proc_ids[i][0]];
743 }
744 }
745
746 delete[] recvCount;
747
748 // Use barrier for error checking
749 Teuchos::reduceAll<int>(comm, Teuchos::REDUCE_MIN, 2, OK, gOK);
750 if (!gOK[0] || !gOK[1]) {
751 delete [] rbuf;
752 delete [] requests;
753 if (!gOK[0])
754 throw std::runtime_error("Max single message length exceeded");
755 else
756 throw std::bad_alloc();
757 }
758
759 gno_t *sbuf = NULL;
760 if (totalsend) sbuf = new gno_t[totalsend];
761 a = 0;
762
763 for(int j = 0; j < num_node_cmaps; j++) {
764 sbuf[a++] = node_cmap_node_cnts[j];
765 for(int i = 0; i < node_cmap_node_cnts[j]; i++) {
766 sbuf[a++] = node_num_map_[node_ids[j][i]-1];
767 sbuf[a++] = sur_elem[node_ids[j][i]-1].size();
768 for(size_t ecnt=0; ecnt < sur_elem[node_ids[j][i]-1].size(); ecnt++) {
769 sbuf[a++] = sur_elem[node_ids[j][i]-1][ecnt];
770 }
771 }
772 }
773
774 delete[] node_cmap_node_cnts;
775 node_cmap_node_cnts = NULL;
776
777 for(int j = 0; j < num_node_cmaps; j++) {
778 delete[] node_ids[j];
779 }
780
781 delete[] node_ids;
782 node_ids = NULL;
783 ArrayRCP<gno_t> sendBuf;
784
785 if (totalsend)
786 sendBuf = ArrayRCP<gno_t>(sbuf, 0, totalsend, true);
787 else
788 sendBuf = Teuchos::null;
789
790 // Send data; can use readySend since receives are posted.
791 offset = 0;
792 for (int i = 0; i < num_node_cmaps; i++) {
793 if (sendCount[node_proc_ids[i][0]]) {
794 try{
795 Teuchos::readySend<int, gno_t>(comm,
796 Teuchos::arrayView(&sendBuf[offset],
797 sendCount[node_proc_ids[i][0]]),
798 node_proc_ids[i][0]);
799 }
801 }
802 offset += sendCount[node_proc_ids[i][0]];
803 }
804
805 for(int j = 0; j < num_node_cmaps; j++) {
806 delete[] node_proc_ids[j];
807 }
808
809 delete[] node_proc_ids;
810 node_proc_ids = NULL;
811 delete[] sendCount;
812
813 // Wait for messages to return.
814 try{
815 Teuchos::waitAll<int>(comm, Teuchos::arrayView(requests, nrecvranks));
816 }
818
819 delete[] requests;
820 a = 0;
821
822 for (int i = 0; i < num_node_cmaps; i++) {
823 int num_nodes_this_processor = rbuf[a++];
824
825 for (int j = 0; j < num_nodes_this_processor; j++) {
826 int this_node = rbuf[a++];
827 int num_elem_this_node = rbuf[a++];
828
829 for (int ncnt = 0; ncnt < num_nodes_; ncnt++) {
830 if (node_num_map_[ncnt] == this_node) {
831 for (int ecnt = 0; ecnt < num_elem_this_node; ecnt++) {
832 sur_elem[ncnt].push_back(rbuf[a++]);
833 }
834
835 break;
836 }
837 }
838 }
839 }
840
841 delete[] rbuf;
842 //}
843
844#ifndef USE_MESH_ADAPTER
845 for(int ecnt=0; ecnt < num_elem_; ecnt++) {
846 eStart_[ecnt] = nEadj_;
847 MapType eAdjMap;
848 int nnodes = nnodes_per_elem;
849 for(int ncnt=0; ncnt < nnodes; ncnt++) {
850 int node = reconnect[ecnt][ncnt]-1;
851 for(size_t i=0; i < sur_elem[node].size(); i++) {
852 int entry = sur_elem[node][i];
853 typename MapType::iterator iter = eAdjMap.find(entry);
854
855 if(element_num_map_[ecnt] != entry &&
856 iter == eAdjMap.end()) {
857 eAdj.push_back(entry);
858 nEadj_++;
859 eAdjMap.insert({entry, entry});
860 }
861 }
862 }
863
864 eAdjMap.clear();
865 }
866#endif
867
868 for(int b = 0; b < num_elem_; b++) {
869 delete[] reconnect[b];
870 }
871
872 delete[] reconnect;
873 reconnect = NULL;
874 eStart_[num_elem_] = nEadj_;
875
876 eAdj_ = new gno_t [nEadj_];
877
878 for (size_t i=0; i < nEadj_; i++) {
879 eAdj_[i] = as<gno_t>(eAdj[i]);
880 }
881
882 delete[] side_nodes;
883 side_nodes = NULL;
884 delete[] mirror_nodes;
885 mirror_nodes = NULL;
886
887 if (nWeightsPerEntity_ > 0) {
888 entityDegreeWeight_ = new bool [nWeightsPerEntity_];
889 for (int i=0; i < nWeightsPerEntity_; i++) {
890 entityDegreeWeight_[i] = false;
891 }
892 }
893}
894
896template <typename User>
898{
899 if (idx >= 0 && idx < nWeightsPerEntity_)
900 entityDegreeWeight_[idx] = true;
901 else
902 std::cout << "WARNING: invalid entity weight index, " << idx << ", ignored"
903 << std::endl;
904}
905
906template <typename User>
908{
909 std::string fn(" PamgenMesh ");
910 std::cout << me << fn
911 << " dim = " << dimension_
912 << " nodes = " << num_nodes_
913 << " nelems = " << num_elem_
914 << std::endl;
915
916 for (int i = 0; i < num_elem_; i++) {
917 std::cout << me << fn << i
918 << " Elem " << element_num_map_[i]
919 << " Coords: ";
920 for (int j = 0; j < dimension_; j++)
921 std::cout << Acoords_[i + j * num_elem_] << " ";
922 std::cout << std::endl;
923 }
924
925#ifndef USE_MESH_ADAPTER
926 for (int i = 0; i < num_elem_; i++) {
927 std::cout << me << fn << i
928 << " Elem " << element_num_map_[i]
929 << " Graph: ";
930 for (int j = eStart_[i]; j < eStart_[i+1]; j++)
931 std::cout << eAdj_[j] << " ";
932 std::cout << std::endl;
933 }
934#endif
935}
936
937} //namespace Zoltan2
938
939#endif
#define Z2_THROW_NOT_IMPLEMENTED
#define Z2_FORWARD_EXCEPTIONS
Forward an exception back through call stack.
Defines the MeshAdapter interface.
This file defines the StridedData class.
MeshAdapter defines the interface for mesh input.
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...
This class represents a mesh.
InputTraits< User >::node_t node_t
PamgenMeshAdapter(const Comm< int > &comm, std::string typestr="region", int nEntWgts=0)
Constructor for mesh with identifiers but no coordinates or edges.
size_t getLocalNumAdjs(MeshEntityType source, MeshEntityType target) const
Returns the number of first 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.
InputTraits< User >::scalar_t scalar_t
size_t getLocalNum2ndAdjs(MeshEntityType sourcetarget, MeshEntityType through) const
if avail2ndAdjs(), returns the number of second adjacencies on this process.
bool areEntityIDsUnique(MeshEntityType etype) const
Provide a pointer to the entity topology types.
size_t getLocalNumOf(MeshEntityType etype) const
Returns the global number of mesh entities of MeshEntityType.
void setWeightIsDegree(int idx)
Specify an index for which the weight should be the degree of the entity \paran idx Zoltan2 will use ...
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 getTopologyViewOf(MeshEntityType etype, enum EntityTopologyType const *&Types) const
Provide a pointer to the entity topology types.
InputTraits< User >::offset_t offset_t
void getCoordinatesViewOf(MeshEntityType etype, const scalar_t *&coords, int &stride, int dim) const
Provide a pointer to one dimension of entity coordinates.
bool availAdjs(MeshEntityType source, MeshEntityType target) const
Returns whether a first adjacency combination is available.
void getIDsViewOf(MeshEntityType etype, const gno_t *&Ids) const
Provide a pointer to this process' identifiers.
void get2ndAdjsView(MeshEntityType sourcetarget, MeshEntityType through, const offset_t *&offsets, const gno_t *&adjacencyIds) const
if avail2ndAdjs(), set pointers to this process' second adjacencies
InputTraits< User >::part_t part_t
bool avail2ndAdjs(MeshEntityType sourcetarget, MeshEntityType through) const
Returns whether a second adjacency combination is available. If combination is not available in the M...
void getAdjsView(MeshEntityType source, MeshEntityType target, const offset_t *&offsets, const gno_t *&adjacencyIds) const
Sets pointers to this process' mesh first adjacencies.
int getDimension() const
Return dimension of the entity coordinates, if any.
map_t::global_ordinal_type gno_t
Created by mbenlioglu on Aug 31, 2020.
EntityTopologyType
Enumerate entity topology types for meshes: points,lines,polygons,triangles,quadrilaterals,...
MeshEntityType
Enumerate entity types for meshes: Regions, Faces, Edges, or Vertices.
static ArrayRCP< ArrayRCP< zscalar_t > > weights
default_offset_t offset_t
The data type to represent offsets.
default_gno_t gno_t
The ordinal type (e.g., int, long, int64_t) that can represent global counts and identifiers.
default_node_t node_t
The Kokkos node type. This is only meaningful for users of Tpetra objects.
default_lno_t lno_t
The ordinal type (e.g., int, long, int64_t) that represents local counts and local indices.
default_part_t part_t
The data type to represent part numbers.
default_scalar_t scalar_t
The data type for weights and coordinates.