Zoltan2
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
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>::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 
386 private:
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 
451 template <typename User>
452 APFMeshAdapter<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 }
718 template <typename User>
719 void 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 
762 template <typename User>
763 void 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
773 template <typename User>
774 void 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 
806 template <typename User>
807 void 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 
855 template <typename User>
856 void 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
virtual void getAdjsView(MeshEntityType source, MeshEntityType target, const offset_t *&offsets, const gno_t *&adjacencyIds) const
Sets pointers to this process&#39; mesh first adjacencies.
InputTraits< User >::scalar_t scalar_t
APFMeshAdapter(const Comm< int > &comm, apf::Mesh *m, std::string primary, std::string adjacency, bool needSecondAdj=false)
virtual bool areEntityIDsUnique(MeshEntityType etype) const
Provide a pointer to the entity topology types.
virtual bool availAdjs(MeshEntityType source, MeshEntityType target) const
Returns whether a first adjacency combination is available.
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 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.
InputTraits< User >::lno_t lno_t
virtual bool avail2ndAdjs(MeshEntityType sourcetarget, MeshEntityType through) const
Returns whether a second adjacency combination is available. If combination is not available in the M...
static ArrayRCP< ArrayRCP< zscalar_t > > weights
virtual void get2ndAdjsView(MeshEntityType sourcetarget, MeshEntityType through, const offset_t *&offsets, const gno_t *&adjacencyIds) const
if avail2ndAdjs(), set pointers to this process&#39; second adjacencies
virtual int getDimension() const
Return dimension of the entity coordinates, if any.
SparseMatrixAdapter_t::part_t part_t
virtual size_t getLocalNum2ndAdjs(MeshEntityType sourcetarget, MeshEntityType through) const
if avail2ndAdjs(), returns the number of second adjacencies on this process.
virtual void getIDsViewOf(MeshEntityType etype, gno_t const *&Ids) const =0
Provide a pointer to this process&#39; identifiers.
virtual void getCoordinatesViewOf(MeshEntityType etype, const scalar_t *&coords, int &stride, int coordDim) const
Provide a pointer to one dimension of entity coordinates.
virtual size_t getLocalNumAdjs(MeshEntityType source, MeshEntityType target) const
Returns the number of first adjacencies on this process.
InputTraits< User >::part_t part_t
EntityTopologyType
Enumerate entity topology types for meshes: points,lines,polygons,triangles,quadrilaterals, polyhedrons, tetrahedrons, hexhedrons, prisms, or pyramids.
enum Zoltan2::EntityTopologyType topologyAPFtoZ2(enum apf::Mesh::Type ttype)
Vector::node_type Node
void applyPartitioningSolution(const User &in, User *&out, const PartitioningSolution< Adapter > &solution) const
Apply a PartitioningSolution to an input.
MeshEntityType
Enumerate entity types for meshes: Regions, Faces, Edges, or Vertices.
virtual int getNumWeightsPerOf(MeshEntityType etype) const
Return the number of weights per entity.
virtual void getTopologyViewOf(MeshEntityType etype, enum EntityTopologyType const *&Types) const
Provide a pointer to the entity topology types.
enum MeshEntityType getPrimaryEntityType() const
Returns the entity to be partitioned, ordered, colored, etc.
InputTraits< User >::offset_t offset_t
This file defines the StridedData class.
Zoltan2::BasicUserTypes< zscalar_t, zlno_t, zgno_t > user_t
Definition: Metric.cpp:74