Zoltan2
Loading...
Searching...
No Matches
Zoltan2_APFMeshAdapter.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_APFMESHADAPTER_HPP_
51#define _ZOLTAN2_APFMESHADAPTER_HPP_
52
55#include <map>
56#include <unordered_map>
57#include <vector>
58#include <string>
59#include <cassert>
60
61#ifndef HAVE_ZOLTAN2_PARMA
62
63namespace apf {
64 class Mesh;
65}
66namespace Zoltan2 {
67template <typename User>
68class APFMeshAdapter : public MeshAdapter<User>
69{
70public:
71
72
73 APFMeshAdapter(const Comm<int> &comm, apf::Mesh* m,std::string primary,std::string adjacency,bool needSecondAdj=false)
74 {
75 throw std::runtime_error(
76 "BUILD ERROR: ParMA requested but not compiled into Zoltan2.\n"
77 "Please set CMake flag Trilinos_ENABLE_SCOREC:BOOL=ON.");
78 }
79};
80}
81#endif
82
83#ifdef HAVE_ZOLTAN2_PARMA
84
85#include <apfMesh.h>
86#include <apfDynamicArray.h>
87#include <apfNumbering.h>
88#include <apfShape.h>
89#include <PCU.h>
90namespace Zoltan2 {
91
118template <typename User>
119 class APFMeshAdapter: public MeshAdapter<User> {
120
121public:
122
123 typedef typename InputTraits<User>::scalar_t scalar_t;
124 typedef typename InputTraits<User>::lno_t lno_t;
125 typedef typename InputTraits<User>::gno_t gno_t;
126 typedef typename InputTraits<User>::offset_t offset_t;
127 typedef typename InputTraits<User>::part_t part_t;
128 typedef typename InputTraits<User>::node_t node_t;
129 typedef User user_t;
130
143 APFMeshAdapter(const Comm<int> &comm, apf::Mesh* m,std::string primary,
144 std::string adjacency,bool needSecondAdj=false, int needs=0);
145
146 void destroy();
147 void print(int me,int verbosity=0);
148 template <typename Adapter>
149 void applyPartitioningSolution(const User &in, User *&out,
150 const PartitioningSolution<Adapter> &solution) const{
151
152 apf::Migration* plan = new apf::Migration(*out);
153 const part_t* new_part_ids = solution.getPartListView();
154
155 if ((m_dimension==3 && this->getPrimaryEntityType()==MESH_REGION) ||
156 (m_dimension==2&&this->getPrimaryEntityType()==MESH_FACE)) {
157 //Elements can simply be sent to the given target parts
158 apf::MeshIterator* itr = (*out)->begin(m_dimension);
159 apf::MeshEntity* ent;
160 int i=0;
161 while ((ent=(*out)->iterate(itr))) {
162 assert(new_part_ids[i]<PCU_Comm_Peers());
163 plan->send(ent,new_part_ids[i]);
164 i++;
165 }
166 }
167 else {
168 //For non-element entities we have to select elements based on the non-element
169 // based Zoltan2 partition. We do this by sending the ith element to the part
170 // that will have the most of the elements downward entities.
171 int dim = entityZ2toAPF(this->getPrimaryEntityType());
172 apf::MeshIterator* itr = (*out)->begin(m_dimension);
173 apf::MeshEntity* ent;
174 size_t i=0;
175 while ((ent=(*out)->iterate(itr))) {
176 std::unordered_map<unsigned int,unsigned int> newOwners;
177 apf::Downward adj;
178 unsigned int max_num = 0;
179 int new_part=PCU_Comm_Self();
180 unsigned int num = in->getDownward(ent,dim,adj);
181 for (unsigned int j=0;j<num;j++) {
182 gno_t gid = apf::getNumber(gids[dim],apf::Node(adj[j],0));
183 lno_t lid = apf::getNumber(lids[dim],adj[j],0,0);
184 newOwners[new_part_ids[lid]]++;
185 if (newOwners[new_part_ids[lid]]>max_num) {
186 max_num=newOwners[new_part_ids[lid]];
187 new_part = new_part_ids[lid];
188 }
189 }
190 if (max_num>1)
191 if (new_part<0||new_part>=PCU_Comm_Peers()) {
192 std::cout<<new_part<<std::endl;
193 throw std::runtime_error("Target part is out of bounds\n");
194 }
195 plan->send(ent,new_part);
196 i++;
197 }
198
199 }
200 (*out)->migrate(plan);
201 }
202
204 // The MeshAdapter interface.
205 // This is the interface that would be called by a model or a problem .
207
208 /* NOTE: Only elements are uniquely provided from the APF Mesh Adapter.
209 All other elements have copies across the shared parts
210 These copies can be joined by the sharing of a unique global id
211 getGlobalNumOf(type) != Sum(getLocalNumOf(type))
212 */
213 bool areEntityIDsUnique(MeshEntityType etype) const {
214 int dim = entityZ2toAPF(etype);
215 return dim==m_dimension;
216 }
217 size_t getLocalNumOf(MeshEntityType etype) const
218 {
219 int dim = entityZ2toAPF(etype);
220 if (dim<=m_dimension&&dim>=0)
221 return num_local[dim];
222 return 0;
223 }
224
225 void getIDsViewOf(MeshEntityType etype, const gno_t *&Ids) const
226 {
227 int dim = entityZ2toAPF(etype);
228 if (dim<=m_dimension&&dim>=0)
229 Ids = gid_mapping[dim];
230 else
231 Ids = NULL;
232 }
233
235 enum EntityTopologyType const *&Types) const {
236 int dim = entityZ2toAPF(etype);
237 if (dim<=m_dimension&&dim>=0)
238 Types = topologies[dim];
239 else
240 Types = NULL;
241 }
242
243 int getNumWeightsPerOf(MeshEntityType etype) const {
244 int dim = entityZ2toAPF(etype);
245 return static_cast<int>(weights[dim].size());
246}
247
248 void getWeightsViewOf(MeshEntityType etype, const scalar_t *&ws,
249 int &stride, int idx = 0) const
250 {
251 int dim = entityZ2toAPF(etype);
252 typename map_array_t::iterator itr = weights[dim].find(idx);
253 if (itr!=weights[dim].end()) {
254 ws = &(*(itr->second.first));
255 stride = itr->second.second;
256 }
257 else {
258 ws = NULL;
259 stride = 0;
260 }
261 }
262
263 int getDimension() const { return coord_dimension; }
264
265 void getCoordinatesViewOf(MeshEntityType etype, const scalar_t *&coords,
266 int &stride, int coordDim) const {
267 if (coordDim>=0 && coordDim<3) {
268 int dim = entityZ2toAPF(etype);
269 if (dim<=m_dimension&&dim>=0) {
270 coords = ent_coords[dim]+coordDim;
271 stride = 3;
272 }
273 else {
274 coords = NULL;
275 stride = 0;
276 }
277 }
278 else {
279 coords = NULL;
280 stride = 0;
281 }
282 }
283
284 bool availAdjs(MeshEntityType source, MeshEntityType target) const {
285 int dim_source = entityZ2toAPF(source);
286 int dim_target = entityZ2toAPF(target);
287 return dim_source<=m_dimension && dim_source>=0 &&
288 dim_target<=m_dimension && dim_target>=0 &&
289 dim_target!=dim_source&&
290 has(dim_source) && has(dim_target);
291 }
292
293 size_t getLocalNumAdjs(MeshEntityType source, MeshEntityType target) const
294 {
295 int dim_source = entityZ2toAPF(source);
296 int dim_target = entityZ2toAPF(target);
297 if (availAdjs(source,target))
298 return adj_gids[dim_source][dim_target].size();
299 return 0;
300 }
301
302 void getAdjsView(MeshEntityType source, MeshEntityType target,
303 const offset_t *&offsets, const gno_t *& adjacencyIds) const
304 {
305 int dim_source = entityZ2toAPF(source);
306 int dim_target = entityZ2toAPF(target);
307 if (availAdjs(source,target)) {
308 offsets = adj_offsets[dim_source][dim_target];
309 adjacencyIds = &(adj_gids[dim_source][dim_target][0]);
310 }
311 else {
312 offsets=NULL;
313 adjacencyIds = NULL;
314 }
315 }
316 //TODO:: some pairings of the second adjacencies do not include off processor adjacencies.
317 // one such pairing is the edge through vertex second adjacnecies.
318 //#define USE_MESH_ADAPTER
319#ifndef USE_MESH_ADAPTER
320 bool avail2ndAdjs(MeshEntityType sourcetarget, MeshEntityType through) const
321 {
322 if (adj2_gids==NULL)
323 return false;
324 int dim_source = entityZ2toAPF(sourcetarget);
325 int dim_target = entityZ2toAPF(through);
326 if (dim_source==1&&dim_target==0)
327 return false;
328 return dim_source<=m_dimension && dim_source>=0 &&
329 dim_target<=m_dimension && dim_target>=0 &&
330 dim_target!=dim_source &&
331 has(dim_source)&&has(dim_target);
332 }
333
334 size_t getLocalNum2ndAdjs(MeshEntityType sourcetarget,
335 MeshEntityType through) const
336 {
337 int dim_source = entityZ2toAPF(sourcetarget);
338 int dim_target = entityZ2toAPF(through);
339 if (avail2ndAdjs(sourcetarget,through))
340 return adj2_gids[dim_source][dim_target].size();
341 return 0;
342
343 }
344
345 void get2ndAdjsView(MeshEntityType sourcetarget, MeshEntityType through,
346 const offset_t *&offsets, const gno_t *&adjacencyIds) const
347 {
348 int dim_source = entityZ2toAPF(sourcetarget);
349 int dim_target = entityZ2toAPF(through);
350 if (avail2ndAdjs(sourcetarget,through)) {
351 offsets=adj2_offsets[dim_source][dim_target];
352 adjacencyIds=&(adj2_gids[dim_source][dim_target][0]);
353 }
354
355 }
356#endif
357
371 void setWeights(MeshEntityType etype, const scalar_t *val, int stride, int idx=0);
372
384 void setWeights(MeshEntityType etype, apf::Mesh* m,apf::MeshTag* weights, int* ids=NULL);
385
386private:
391 bool has(int dim) const {return (entity_needs>>dim)%2;}
392
393 // provides a conversion from the mesh entity type to the apf dimension
394 int entityZ2toAPF(enum MeshEntityType etype) const {return static_cast<int>(etype);}
395
396 // provides a conversion from the apf topology type to the Zoltan2 topology type
397 enum EntityTopologyType topologyAPFtoZ2(enum apf::Mesh::Type ttype) const {
398 if (ttype==apf::Mesh::VERTEX)
399 return POINT;
400 else if (ttype==apf::Mesh::EDGE)
401 return LINE_SEGMENT;
402 else if (ttype==apf::Mesh::TRIANGLE)
403 return TRIANGLE;
404 else if (ttype==apf::Mesh::QUAD)
405 return QUADRILATERAL;
406 else if (ttype==apf::Mesh::TET)
407 return TETRAHEDRON;
408 else if (ttype==apf::Mesh::HEX)
409 return HEXAHEDRON;
410 else if (ttype==apf::Mesh::PRISM)
411 return PRISM;
412 else if (ttype==apf::Mesh::PYRAMID)
413 return PYRAMID;
414 else
415 throw "No such APF topology type";
416
417 }
418
419 // provides a conversion from the mesh tag type to scalar_t since mesh tags are not templated
420 void getTagWeight(apf::Mesh* m, apf::MeshTag* tag,apf::MeshEntity* ent, scalar_t* ws);
421
422
423 int m_dimension; //Dimension of the mesh
424
425 //An int between 0 and 15 that represents the mesh dimensions that are constructed
426 // in binary. A 1 in the ith digit corresponds to the ith dimension being constructed
427 // Ex: 9 = 1001 is equivalent to regions and vertices are needed
428 int entity_needs;
429 apf::Numbering** lids; //[dimension] numbering of local id numbers
430 apf::GlobalNumbering** gids;//[dimension] numbering of global id numbers
431 gno_t** gid_mapping; //[dimension][lid] corresponding global id numbers
432 size_t* num_local; //[dimension] number of local entities
433 EntityTopologyType** topologies; //[dimension] topologies for each entity
434 offset_t*** adj_offsets; //[first_dimension][second_dimension] array of offsets
435 std::vector<gno_t>** adj_gids; //[first_dimension][second_dimension] global_ids of first adjacencies
436 offset_t*** adj2_offsets; //[first_dimension][second_dimension] array of offsets for second adjacencies
437 std::vector<gno_t>** adj2_gids; //[first_dimension][second_dimension] global_ids of second adjacencies
438 int coord_dimension; //dimension of coordinates (always 3 for APF)
439 scalar_t** ent_coords; //[dimension] array of coordinates [xs ys zs]
440
441 //[dimension][id] has the start of the weights array and the stride
442 typedef std::unordered_map<int, std::pair<ArrayRCP<const scalar_t>, int> > map_array_t;
443 map_array_t* weights;
444
445};
446
448// Definitions
450
451template <typename User>
452APFMeshAdapter<User>::APFMeshAdapter(const Comm<int> &comm,
453 apf::Mesh* m,
454 std::string primary,
455 std::string adjacency,
456 bool needSecondAdj,
457 int needs) {
458
459 //get the mesh dimension
460 m_dimension = m->getDimension();
461
462 //get the dimensions that are needed to be constructed
463 entity_needs = needs;
464
465 //Make the primary and adjacency entity types
466 //choices are region, face, edge, vertex
467 //element is a shortcut to mean the mesh dimension entity type
468 //region will throw an error on 2D meshes
469 if (primary=="element") {
470 if (m_dimension==2)
471 primary="face";
472 else
473 primary="region";
474 }
475 if (adjacency=="element") {
476 if (m_dimension==2)
477 adjacency="face";
478 else
479 adjacency="region";
480 }
481 if (primary=="region"&&m_dimension<3)
482 throw std::runtime_error("primary type and mesh dimension mismatch");
483 if (adjacency=="region"&&m_dimension<3)
484 throw std::runtime_error("adjacency type and mesh dimension mismatch");
485 this->setEntityTypes(primary,adjacency,adjacency);
486
487 //setup default needs such that primary and adjacency types are always constructed
488 int dim1 = entityZ2toAPF(this->getPrimaryEntityType());
489 int dim2 = entityZ2toAPF(this->getAdjacencyEntityType());
490 int new_needs=0;
491 new_needs+=1<<dim1;
492 new_needs+=1<<dim2;
493 entity_needs|=new_needs;
494
495 //count the local and global numbers as well as assign ids and map local to global
496 lids = new apf::Numbering*[m_dimension+1];
497 gids = new apf::GlobalNumbering*[m_dimension+1];
498 gid_mapping = new gno_t*[m_dimension+1];
499 std::unordered_map<gno_t,lno_t>* lid_mapping = new std::unordered_map<gno_t,lno_t>[m_dimension+1];
500 num_local = new size_t[m_dimension+1];
501 topologies = new EntityTopologyType*[m_dimension+1];
502
503 for (int i=0;i<=m_dimension;i++) {
504 num_local[i]=0;
505
506 topologies[i] = NULL;
507 gid_mapping[i] = NULL;
508 if (!has(i))
509 continue;
510 //number of local and global entities
511 num_local[i] = m->count(i);
512 long global_count = countOwned(m,i);
513 PCU_Add_Longs(&global_count,1);
514
515
516 //Number each entity with local and global numbers
517 char lids_name[15];
518 sprintf(lids_name,"lids%d",i);
519 char gids_name[15];
520 sprintf(gids_name,"ids%d",i);
521 apf::FieldShape* shape = apf::getConstant(i);
522 lids[i] = apf::createNumbering(m,lids_name,shape,1);
523 apf::Numbering* tmp = apf::numberOwnedDimension(m,gids_name,i);
524 gids[i] = apf::makeGlobal(tmp);
525 apf::synchronize(gids[i]);
526 apf::MeshIterator* itr = m->begin(i);
527 apf::MeshEntity* ent;
528 unsigned int num=0;
529 while ((ent=m->iterate(itr))) {
530 apf::number(lids[i],ent,0,0,num);
531 lid_mapping[i][apf::getNumber(gids[i],apf::Node(ent,0))]=num;
532 num++;
533 }
534 m->end(itr);
535 assert(num==num_local[i]);
536
537 //Make a mapping from local to global
538 //While we are at it take the topology types
539 gid_mapping[i] = new gno_t[num_local[i]];
540 topologies[i] = new EntityTopologyType[num_local[i]];
541 apf::DynamicArray<apf::Node> nodes;
542 itr = m->begin(i);
543 num=0;
544 while((ent=m->iterate(itr))) {
545 gno_t gid = apf::getNumber(gids[i],apf::Node(ent,0));
546 gid_mapping[i][ apf::getNumber(lids[i],ent,0,0)] = gid;
547 topologies[i][num] = topologyAPFtoZ2(m->getType(ent));
548 num++;
549 }
550 m->end(itr);
551
552
553 }
554 //First Adjacency and Second Adjacency data
555 adj_gids = new std::vector<gno_t>*[m_dimension+1];
556 adj_offsets = new offset_t**[m_dimension+1];
557 if (needSecondAdj) {
558 adj2_gids = new std::vector<gno_t>*[m_dimension+1];
559 adj2_offsets = new offset_t**[m_dimension+1];
560 }
561 else {
562 adj2_gids=NULL;
563 adj2_offsets=NULL;
564 }
565 for (int i=0;i<=m_dimension;i++) {
566 adj_gids[i]=NULL;
567 adj_offsets[i]=NULL;
568 if (needSecondAdj) {
569 adj2_gids[i]=NULL;
570 adj2_offsets[i]=NULL;
571 }
572 if (!has(i))
573 continue;
574 adj_gids[i] = new std::vector<gno_t>[m_dimension+1];
575 adj_offsets[i] = new offset_t*[m_dimension+1];
576 if (needSecondAdj) {
577 adj2_gids[i] = new std::vector<gno_t>[m_dimension+1];
578 adj2_offsets[i] = new offset_t*[m_dimension+1];
579 }
580 for (int j=0;j<=m_dimension;j++) {
581
582 if (i==j||!has(j)) {
583 adj_offsets[i][j]=NULL;
584 if (needSecondAdj)
585 adj2_offsets[i][j]=NULL;
586 continue;
587 }
588
589 //Loop through each entity
590 apf::MeshIterator* itr = m->begin(i);
591 apf::MeshEntity* ent;
592
593 adj_offsets[i][j] = new offset_t[num_local[i]+1];
594 adj_offsets[i][j][0] =0;
595 if (needSecondAdj) {
596 adj2_offsets[i][j] = new offset_t[num_local[i]+1];
597 adj2_offsets[i][j][0] =0;
598 }
599 int k=1;
600
601 //We need communication for second adjacency
602 if (needSecondAdj)
603 PCU_Comm_Begin();
604 std::unordered_map<gno_t,apf::MeshEntity*> part_boundary_mapping;
605
606 while ((ent=m->iterate(itr))) {
607 std::set<gno_t> temp_adjs; //temp storage for second adjacency
608 //Get First Adjacency
609 apf::Adjacent adj;
610 m->getAdjacent(ent,j,adj);
611 for (unsigned int l=0;l<adj.getSize();l++) {
612 adj_gids[i][j].push_back(apf::getNumber(gids[j],apf::Node(adj[l],0)));
613 //Now look at Second Adjacency
614 if (needSecondAdj) {
615 apf::Adjacent adj2;
616 m->getAdjacent(adj[l],i,adj2);
617 for (unsigned int o=0;o<adj2.getSize();o++)
618 temp_adjs.insert(apf::getNumber(gids[i],apf::Node(adj2[o],0)));
619 if (i==m_dimension) {
620 apf::Parts res;
621 m->getResidence(adj[l],res);
622
623 part_boundary_mapping[apf::getNumber(gids[j],apf::Node(adj[l],0))] = adj[l];
624 for (apf::Parts::iterator it=res.begin();it!=res.end();it++) {
625 gno_t send_vals[2];
626 send_vals[1]=apf::getNumber(gids[i],apf::Node(ent,0));
627 send_vals[0]=apf::getNumber(gids[j],apf::Node(adj[l],0));
628
629 PCU_Comm_Pack(*it,send_vals,2*sizeof(gno_t));
630 }
631 }
632 }
633 }
634 adj_offsets[i][j][k] = adj_gids[i][j].size();
635 k++;
636 //Copy over local second adjacencies to copies
637 if (needSecondAdj && i!=m_dimension) {
638 apf::Parts res;
639 m->getResidence(ent,res);
640 typename std::set<gno_t>::iterator adj_itr;
641 for (adj_itr=temp_adjs.begin();adj_itr!=temp_adjs.end();adj_itr++) {
642 for (apf::Parts::iterator it=res.begin();it!=res.end();it++) {
643 gno_t send_vals[2];
644 send_vals[0]=apf::getNumber(gids[i],apf::Node(ent,0));
645 send_vals[1] = *adj_itr;
646 if (send_vals[0]!=send_vals[1])
647 PCU_Comm_Pack(*it,send_vals,2*sizeof(gno_t));
648 }
649 }
650 }
651 }
652 m->end(itr);
653 if (needSecondAdj) {
654 //Now capture mesh wide second adjacency locally
655 PCU_Comm_Send();
656 std::set<gno_t>* adjs2 = new std::set<gno_t>[num_local[i]];
657 while (PCU_Comm_Receive()) {
658 gno_t adj2[2];
659 PCU_Comm_Unpack(adj2,2*sizeof(gno_t));
660
661 if (i==m_dimension) {
662 apf::MeshEntity* through = part_boundary_mapping[adj2[0]];
663 apf::Adjacent adj;
664 m->getAdjacent(through,i,adj);
665 for (unsigned int l=0;l<adj.getSize();l++) {
666 if (apf::getNumber(gids[i],apf::Node(adj[l],0))!=adj2[1])
667 adjs2[apf::getNumber(lids[i],adj[l],0,0)].insert(adj2[1]);
668 }
669 }
670 else {
671 lno_t index = lid_mapping[i][adj2[0]];
672 adjs2[index].insert(adj2[1]);
673 }
674 }
675 //And finally convert the second adjacency to a vector to be returned to user
676 for (size_t l=0;l<num_local[i];l++) {
677 for (typename std::set<gno_t>::iterator sitr = adjs2[l].begin();sitr!=adjs2[l].end();sitr++) {
678 adj2_gids[i][j].push_back(*sitr);
679 }
680 adj2_offsets[i][j][l+1]=adj2_gids[i][j].size();
681 }
682 }
683 }
684 }
685 //Coordinates
686 coord_dimension = 3;
687 ent_coords = new scalar_t*[m_dimension+1];
688 for (int i=0;i<=m_dimension;i++) {
689 ent_coords[i] = NULL;
690 if (!has(i))
691 continue;
692 apf::MeshIterator* itr = m->begin(i);
693 apf::MeshEntity* ent;
694 ent_coords[i] = new scalar_t[3*num_local[i]];
695 int j=0;
696 while((ent=m->iterate(itr))) {
697 apf::Vector3 point;
698 if (i==0) {
699 m->getPoint(ent,0,point);
700 }
701 else {
702 point = apf::getLinearCentroid(m,ent);
703 }
704 for (int k=0;k<3;k++)
705 ent_coords[i][j*3+k] = point[k];
706 j++;
707 }
708 m->end(itr);
709 }
710
711 //Just make the weights array with nothing in it for now
712 //It will be filled by calls to setWeights(...)
713 weights = new map_array_t[m_dimension+1];
714
715 //cleanup
716 delete [] lid_mapping;
717}
718template <typename User>
719void APFMeshAdapter<User>::destroy() {
720 //So that we can't destory the adapter twice
721 if (m_dimension==-1)
722 return;
723 for (int i=0;i<=m_dimension;i++) {
724 if (!has(i))
725 continue;
726 delete [] ent_coords[i];
727 delete [] adj_gids[i];
728 if (adj2_gids)
729 delete [] adj2_gids[i];
730 for (int j=0;j<=m_dimension;j++) {
731 if (!has(j))
732 continue;
733 if (i!=j) {
734 delete [] adj_offsets[i][j];
735 if (adj2_gids)
736 delete [] adj2_offsets[i][j];
737 }
738 }
739 if (adj2_gids)
740 delete [] adj2_offsets[i];
741 delete [] adj_offsets[i];
742 delete [] gid_mapping[i];
743 apf::destroyGlobalNumbering(gids[i]);
744 apf::destroyNumbering(lids[i]);
745 }
746 delete [] ent_coords;
747 delete [] adj_gids;
748 delete [] adj_offsets;
749 if (adj2_gids) {
750 delete [] adj2_gids;
751 delete [] adj2_offsets;
752 }
753 delete [] gid_mapping;
754 delete [] gids;
755 delete [] lids;
756 delete [] num_local;
757 delete [] weights;
758 //Set the mesh dimension to -1 so that no operations can be done on the destroyed adapter
759 m_dimension=-1;
760}
761
762template <typename User>
763void APFMeshAdapter<User>::setWeights(MeshEntityType etype, const scalar_t *val, int stride, int idx) {
764 int dim = entityZ2toAPF(etype);
765 if (dim>m_dimension||!has(dim)) {
766 throw std::runtime_error("Cannot add weights to non existing dimension");
767 }
768 ArrayRCP<const scalar_t> weight_rcp(val,0,stride*getLocalNumOf(etype),false);
769 weights[dim][idx] =std::make_pair(weight_rcp,stride);
770}
771
772//Simple helper function to convert the tag type to the scalar_t type
773template <typename User>
774void APFMeshAdapter<User>::getTagWeight(apf::Mesh* m,
775 apf::MeshTag* tag,
776 apf::MeshEntity* ent,
777 scalar_t* ws) {
778 int size = m->getTagSize(tag);
779 int type = m->getTagType(tag);
780 if (type==apf::Mesh::DOUBLE) {
781 double* w = new double[size];
782 m->getDoubleTag(ent,tag,w);
783 for (int i=0;i<size;i++)
784 ws[i] = static_cast<scalar_t>(w[i]);
785 delete [] w;
786 }
787 else if (type==apf::Mesh::INT) {
788 int* w = new int[size];
789 m->getIntTag(ent,tag,w);
790 for (int i=0;i<size;i++)
791 ws[i] = static_cast<scalar_t>(w[i]);
792 delete [] w;
793 }
794 else if (type==apf::Mesh::LONG) {
795 long* w = new long[size];
796 m->getLongTag(ent,tag,w);
797 for (int i=0;i<size;i++)
798 ws[i] = static_cast<scalar_t>(w[i]);
799 delete [] w;
800 }
801 else {
802 throw std::runtime_error("Unrecognized tag type");
803 }
804}
805
806template <typename User>
807void APFMeshAdapter<User>::setWeights(MeshEntityType etype, apf::Mesh* m,apf::MeshTag* tag, int* ids) {
808 int dim = entityZ2toAPF(etype);
809 if (dim>m_dimension||!has(dim)) {
810 throw std::runtime_error("Cannot add weights to non existing dimension");
811 }
812 int n_weights = m->getTagSize(tag);
813 bool delete_ids = false;
814 if (ids==NULL) {
815 ids = new int[n_weights];
816 delete_ids=true;
817 for (int i=0;i<n_weights;i++)
818 ids[i] = i;
819 }
820 scalar_t* ones = new scalar_t[n_weights];
821 for (int i=0;i<n_weights;i++)
822 ones[i] = 1;
823
824 scalar_t* ws = new scalar_t[num_local[dim]*n_weights];
825 apf::MeshIterator* itr = m->begin(dim);
826 apf::MeshEntity* ent;
827 int j=0;
828 while ((ent=m->iterate(itr))) {
829 scalar_t* w;
830 if (m->hasTag(ent,tag)) {
831 w = new scalar_t[n_weights];
832 getTagWeight(m,tag,ent,w);
833 }
834 else
835 w = ones;
836
837 for (int i=0;i<n_weights;i++) {
838 ws[i*getLocalNumOf(etype)+j] = w[i];
839 }
840 j++;
841
842 if (m->hasTag(ent,tag))
843 delete [] w;
844 }
845 for (int i=0;i<n_weights;i++) {
846 ArrayRCP<const scalar_t> weight_rcp(ws+i*getLocalNumOf(etype),0,getLocalNumOf(etype),i==0);
847 weights[dim][ids[i]] =std::make_pair(weight_rcp,1);
848 }
849
850 if (delete_ids)
851 delete [] ids;
852 delete [] ones;
853}
854
855template <typename User>
856void APFMeshAdapter<User>::print(int me,int verbosity)
857{
858 if (m_dimension==-1) {
859 std::cout<<"Cannot print destroyed mesh adapter\n";
860 return;
861 }
862
863 std::string fn(" APFMesh ");
864 std::cout << me << fn
865 << " dimension = " << m_dimension
866 << std::endl;
867 if (verbosity==0)
868 return;
869 for (int i=0;i<=m_dimension;i++) {
870 if (!has(i))
871 continue;
872 std::cout<<me<<" Number of dimension " << i<< " = " <<num_local[i] <<std::endl;
873 if (verbosity>=1) {
874 for (size_t j=0;j<num_local[i];j++) {
875 std::cout<<" Entity "<<gid_mapping[i][j]<<"("<<j<<"):\n";
876 for (int k=0;k<=m_dimension;k++) {
877 if (!has(k))
878 continue;
879 if (k==i)
880 continue;
881 std::cout<<" First Adjacency of Dimension "<<k<<":";
882 for (offset_t l=adj_offsets[i][k][j];l<adj_offsets[i][k][j+1];l++)
883 std::cout<<" "<<adj_gids[i][k][l];
884 std::cout<<"\n";
885 if (verbosity>=3) {
886 std::cout<<" Second Adjacency through Dimension "<<k<<":";
887 for (offset_t l=adj2_offsets[i][k][j];l<adj2_offsets[i][k][j+1];l++)
888 std::cout<<" "<<adj2_gids[i][k][l];
889 std::cout<<"\n";
890 }
891 }
892 }
893 }
894 }
895}
896
897
898
899} //namespace Zoltan2
900
901#endif //HAVE_ZOLTAN2_PARMA
902
903#endif
enum Zoltan2::EntityTopologyType topologyAPFtoZ2(enum apf::Mesh::Type ttype)
Zoltan2::BasicUserTypes< zscalar_t, zlno_t, zgno_t > user_t
Definition Metric.cpp:74
Defines the MeshAdapter interface.
This file defines the StridedData class.
APFMeshAdapter(const Comm< int > &comm, apf::Mesh *m, std::string primary, std::string adjacency, bool needSecondAdj=false)
InputTraits< User >::node_t node_t
InputTraits< User >::offset_t offset_t
InputTraits< User >::part_t part_t
InputTraits< User >::scalar_t scalar_t
InputTraits< User >::lno_t lno_t
void applyPartitioningSolution(const User &in, User *&out, const PartitioningSolution< Adapter > &solution) const
Apply a PartitioningSolution to an input.
InputTraits< User >::gno_t gno_t
MeshAdapter defines the interface for mesh input.
virtual 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.
virtual bool availAdjs(MeshEntityType source, MeshEntityType target) const
Returns whether a first adjacency combination is available.
virtual int getNumWeightsPerOf(MeshEntityType etype) const
Return the number of weights per entity.
virtual size_t getLocalNum2ndAdjs(MeshEntityType sourcetarget, MeshEntityType through) const
if avail2ndAdjs(), returns the number of second adjacencies on this process.
virtual void getTopologyViewOf(MeshEntityType etype, enum EntityTopologyType const *&Types) const
Provide a pointer to the entity topology types.
virtual int getDimension() const
Return dimension of the entity coordinates, if any.
virtual size_t getLocalNumOf(MeshEntityType etype) const =0
Returns the global number of mesh entities of MeshEntityType.
virtual size_t getLocalNumAdjs(MeshEntityType source, MeshEntityType target) const
Returns the number of first adjacencies on this process.
virtual void getIDsViewOf(MeshEntityType etype, gno_t const *&Ids) const =0
Provide a pointer to this process' identifiers.
virtual bool areEntityIDsUnique(MeshEntityType etype) const
Provide a pointer to the entity topology types.
virtual void getCoordinatesViewOf(MeshEntityType etype, const scalar_t *&coords, int &stride, int coordDim) const
Provide a pointer to one dimension of entity coordinates.
virtual void getAdjsView(MeshEntityType source, MeshEntityType target, const offset_t *&offsets, const gno_t *&adjacencyIds) const
Sets pointers to this process' mesh first adjacencies.
virtual bool avail2ndAdjs(MeshEntityType sourcetarget, MeshEntityType through) const
Returns whether a second adjacency combination is available. If combination is not available in the M...
enum MeshEntityType getPrimaryEntityType() const
Returns the entity to be partitioned, ordered, colored, etc.
virtual void get2ndAdjsView(MeshEntityType sourcetarget, MeshEntityType through, const offset_t *&offsets, const gno_t *&adjacencyIds) const
if avail2ndAdjs(), set pointers to this process' second adjacencies
map_t::local_ordinal_type lno_t
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
SparseMatrixAdapter_t::part_t part_t
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.