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