OpenFPM_pdata  4.1.0 Project that contain the implementation of distributed structures
gargabe.hpp
1 /*
2  * gargabe.hpp
3  *
4  * Created on: Jan 13, 2015
5  * Author: i-bird
6  */
7
8 #ifdef GARGABE_HPP_
9 #define GARGABE_HPP_
10
11
12
13  template <unsigned int j, unsigned int i, typename Graph> void optimize(size_t start_p, Graph & graph)
14  {
15  // We assume that Graph is the rapresentation of a cartesian graph
16  // this mean that the direction d is at the child d
17
18  // Create an Hyper-cube
19
20  HyperCube<dim> hyp;
21
22  // Get the number of wavefronts
23
24  size_t n_wf = hyp.getNumberOfElements_R(0);
25
26  // Get the number of intersecting wavefront
27
28
29
30  // Get the number of sub-dimensional common wavefront
31  // basically are a list of all the subdomain common to two or more
32
33  // Create n_wf wavefront queue
34
36  v.reserve(n_wf);
37
38  // direction of expansion
39
40  size_t domain_id = 0;
41  int exp_dir = 0;
42  bool can_expand = true;
43
44  // while is possible to expand
45
46  while (can_expand)
47  {
48  // for each direction of expansion expand the wavefront
49
50  for (int d = 0 ; d < n_wf ; d++)
51  {
52  // get the wavefront at direction d
53
54  openfpm::vector<size_t> & wf_d = v_w.get<wavefront::domains>(d);
55
56  // flag to indicate if the wavefront can expand
57
58  bool w_can_expand = true;
59
60  // for each subdomain
61
62  for (size_t sub = 0 ; sub < wf_d.size() ; sub++)
63  {
64  // check if the adjacent domain in direction d exist
65  // and is of the same id
66
67  // get the starting subdomain
68  size_t sub_w = wf_d.get<0>(sub);
69
70  // we get the processor id of the neighborhood sub-domain on direction d
71  size_t exp_p = graph.getChild(sub_w,d).get<j>();
72
73  // we check if it is the same processor id
74  if (exp_p != domain_id)
75  {
76  w_can_expand = false;
77  }
78  }
79
80  // if we can expand the wavefront expand it
81  if (w_can_expand == true)
82  {
83  // for each subdomain
84  for (size_t sub = 0 ; sub < wf_d.size() ; sub++)
85  {
86  // update the position of the wavefront
87  wf_d.get<0>(sub) = wf_d.get<0>(sub) + gh.stride(d);
88  }
89
90  // here we add sub-domains to all the other queues
91  // get the face of the hyper-cube
92
93  SubHyperCube<dim,dim-1> sub_hyp = hyp.getSubHyperCube(d);
94
95  std::vector<comb<dim>> q_comb = sub_hyp.getCombinations_R(dim-2);
96  }
97  }
98  }
99
100  // For each point in the Hyper-cube check if we can move the wave front
101
102
103  }
104
105 #ifndef PARALLEL_DECOMPOSITION
106 // CreateSubspaces();
107 #endif
108
109 #ifndef USE_METIS_GP
110
111  // Here we do not use METIS
112  // Distribute the divided domains
113
114  // Get the number of processing units
115  size_t Np = v_cl.getProcessingUnits();
116
117  // Get the ID of this processing unit
118  // and push the subspace is taking this
119  // processing unit
120
121  for (size_t p_id = v_cl.getProcessUnitID(); p_id < Np ; p_id += Np)
122  id_sub.push_back(p_id);
123 #else
124
125
126 #endif
127
128
129
132
133  // get the decomposition
134  auto & dec = g_dist.getDecomposition();
135
136  Vcluster & v_cl = *global_v_cluster;
137
138  // check the consistency of the decomposition
139  val = dec.check_consistency();
140  BOOST_REQUIRE_EQUAL(val,true);
141
142  // for each local volume
143  // Get the number of local grid needed
144  size_t n_grid = dec.getNLocalHyperCube();
145
146  size_t vol = 0;
147
149
150  // Allocate the grids
151  for (size_t i = 0 ; i < n_grid ; i++)
152  {
153  // Get the local hyper-cube
154  SpaceBox<2,float> sub = dec.getLocalHyperCube(i);
155
156  Box<2,size_t> g_box = g_dist.getCellDecomposer().convertDomainSpaceIntoGridUnits(sub);
158
159  vol += g_box.getVolumeKey();
160  }
161
162  v_cl.reduce(vol);
163  v_cl.execute();
164
165  BOOST_REQUIRE_EQUAL(vol,k*k);
166
168
169
170  // 3D test
171
172  // g_dist.write("");
173
174  /* auto g_it = g_dist.getIteratorBulk();
175
176  auto g_it_halo = g_dist.getHalo();
177
178  // Let try to solve the poisson equation d2(u) = f with f = 1 and computation
179  // comunication overlap (100 Jacobi iteration)
180
181  for (int i = 0 ; i < 100 ; i++)
182  {
183  g_dist.ghost_get();
184
185  // Compute the bulk
186
187  jacobi_iteration(g_it);
188
189  g_dist.ghost_sync();
190
191  // Compute the halo
192
193  jacobi_iteration(g_it_halo);
194  }*/
195
196
197  BOOST_AUTO_TEST_CASE( grid_dist_id_poisson_test_use)
198  {
199  // grid size
200  /* size_t sz[2] = {1024,1024};
201
202  // Distributed grid with id decomposition
203
204  grid_dist_id<2, scalar<float>, CartDecomposition<2,size_t>> g_dist(sz);
205
206  // Create the grid on memory
207
208  g_dist.Create();*/
209
210  /* auto g_it = g_dist.getIteratorBulk();
211
212  auto g_it_halo = g_dist.getHalo();
213
214  // Let try to solve the poisson equation d2(u) = f with f = 1 and computation
215  // comunication overlap (100 Jacobi iteration)
216
217  for (int i = 0 ; i < 100 ; i++)
218  {
219  g_dist.ghost_get();
220
221  // Compute the bulk
222
223  jacobi_iteration(g_it);
224
225  g_dist.ghost_sync();
226
227  // Compute the halo
228
229  jacobi_iteration(g_it_halo);
230  }*/
231  }
232
233  template<typename iterator> void jacobi_iteration(iterator g_it, grid_dist_id<2, float, scalar<float>, CartDecomposition<2,float>> & g_dist)
234  {
235  // scalar
236  typedef scalar<float> S;
237
238  // iterator
239
240  while(g_it.isNext())
241  {
242  // Jacobi update
243
244  auto pos = g_it.get();
245
246  g_dist.template get<S::ele>(pos) = (g_dist.template get<S::ele>(pos.move(0,1)) +
247  g_dist.template get<S::ele>(pos.move(0,-1)) +
248  g_dist.template get<S::ele>(pos.move(1,1)) +
249  g_dist.template get<S::ele>(pos.move(1,-1)) / 4.0);
250
251  ++g_it;
252  }
253  }
254
255 =======
256
257  /*
258  * CartDecomposition.cpp
259  *
260  * Created on: Aug 15, 2014
261  * Author: Pietro Incardona
262  */
263
264  #include "CartDecomposition.hpp"
265
266
267
276  /*template<typename T> T CartDecomposition<T>::getBulk(T data)
277  {
278  // for each element in data
279
280  for (size_t i = 0; i < data.size() ; i++)
281  {
282  if (localSpace.isInside())
283  }
284
285  }
286
287  template<typename T> T CartDecomposition<T>::getInternal()
288  {
289
290  }*/
291
299  bool borderOrBulk(neighborhood & nb)
300  {
301  device::grid<1,size_t> nbr = nb.next();
302
303  // check the neighborhood
304
305  // get neighborhood iterator
306
307  grid_key_dx_iterator<dim> iterator_nbr = nbr.getIterator();
308
309  while (iterator_nbr.hasNext())
310  {
311  grid_key_dx key_nbr = iterator_nbr.next();
312
313  // check if the neighboorhood is internal
314
315  if(subspace.isBound(data.template get<Point::x>(key_nbr)) == false)
316  {
317  // it is border
318
319  return true;
320
321  ret.bord.push_back(key);
322  break;
323  }
324  }
325
326  return false;
327  }
328
341  template<unsigned int dim, typename T, template<typename> class layout, typename Memory, template<unsigned int, typename> class Domain, template<typename, typename, typename> class data_s>
342  dataDiv<T> CartDecomposition<dim,T,layout>::divide(device::grid<1,Point<dim,T>> & data, neighborhood & nb)
343  {
345
347
348  ret.bord = new boost::shared_ptr<T>(new T());
349  ret.inte = new boost::shared_ptr<T>(new T());
350  ret.ext = new boost::shared_ptr<T>(new T());
351
353
354  grid_key_dx_iterator<dim> iterator = data.getIterator();
355
357
358  while (iterator.hasNext())
359  {
360  grid_key_dx<dim> key = iterator.next();
361
363
364  if (subspace.isBound(data.template get<Point<3,T>::x>(key)))
365  {
367
368  if (borderOrBulk(nb) == true)
369  {
370  // It is border
371
372  ret.bord.push_back(key);
373  }
374  else
375  {
376  // It is bulk
377
378  ret.bulk.push_back(key);
379  }
380  }
381  else
382  {
384
385  ret.ext.push_back(key);
386  }
387  }
388  }
389
390
391 >>>>>>> Jenkin script for taurus
392
393
402 /* template <typename obj> Vcluster_object_array<obj> allocate(size_t n)
403 {
404  // Vcluster object array
405  Vcluster_object_array<obj> vo;
406
407  // resize the array
408  vo.resize(n);
409
410  // Create the object on memory and return a Vcluster_object_array
411  return vo;
412 }*/
413
414
415 /*template<typename T>
416 class Vcluster_object_array : public VObject
417 {
418  std::vector<T> objects;
419
420 public:*/
421
425 /* Vcluster_object_array()
426  {
427
428  }*/
429
435 /* size_t size() const
436  {
437  return objects.size();
438  }*/
439
446 /* T & get(unsigned int i)
447  {
448  return objects[i];
449  }*/
450
456 /* const T & get(unsigned int i) const
457  {
458  return objects[i];
459  }*/
460
466 /* bool isArray()
467  {
468  return true;
469  }*/
470
474 /* virtual void destroy()
475  {
476  // Destroy the objects
477  objects.clear();
478  }*/
479
485 /* size_t packObjectSize()
486  {
487  size_t message = 0;
488
489  // Destroy each objects
490  for (size_t i = 0 ; i < objects.size() ; i++)
491  {
492  message += objects[i].packObjectSize();
493  }
494
495  return message;
496  }*/
497
498
506 /* size_t packObject(void * mem)
507  {
508  // Pointer is zero
509  size_t ptr = 0;
510  unsigned char * m = (unsigned char *)mem;
511
512  // pack each object
513  for (size_t i = 0 ; i < objects.size() ; i++)
514  {
515  ptr += objects[i].packObject(&m[ptr]);
516  }
517
518 #ifdef DEBUG
519  if (ptr != packObjectSize())
520  {
521  std::cerr << "Error " << __FILE__ << " " << __LINE__ << " the pack object size does not match the message" << "\n";
522  }
523 #endif
524
525  return ptr;
526  }*/
527
533 /* size_t packObjectInArraySize(size_t i)
534  {
535  return objects[i].packObjectSize();
536  }*/
537
545 /* size_t packObjectInArray(size_t i, void * p)
546  {
547  return objects[i].packObject(p);
548  }*/
549
555 /* void destroy(size_t i)
556  {
557  objects.erase(objects.begin() + i);
558  }*/
559
565 /* T & operator[](size_t j)
566  {
567  return objects[j];
568  }*/
569
575 /* const T & operator[](size_t j) const
576  {
577  return objects[j];
578  }*/
579
585 /* void resize(size_t n)
586  {
587  objects.resize(n);
588  }
589 };*/
590
597 /*class VObject
598 {
599 public:
600
601  // Check if this Object is an array
602  virtual bool isArray() = 0;
603
604  // destroy the object
605  virtual void destroy() = 0;
606
607  // get the size of the memory needed to pack the object
608  virtual size_t packObjectSize() = 0;
609
610  // pack the object
611  virtual size_t packObject(void *) = 0;
612
613  // get the size of the memory needed to pack the object in the array
614  virtual size_t packObjectInArraySize(size_t i) = 0;
615
616  // pack the object in the array (the message produced can be used to move one)
617  // object from one processor to another
618  virtual size_t packObjectInArray(size_t i, void * p) = 0;
619
620  // destroy an element from the array
621  virtual void destroy(size_t n) = 0;
622 };*/
623
624
625
626
627
628
629
630
631
646 template<typename T> void impose(const T & op , typename Sys_eqs::stype num ,long int id ,grid_dist_iterator_sub<Sys_eqs::dims,typename g_map_type::d_grid> it_d, bool skip_first = false)
647 {
649
652
653  // Construct the map 3 processors 1 processors
654
655  std::unordered_map<size_t,size_t> map_row;
656
657  auto it2 = g_map.getDomainGhostIterator();
658  auto ginfo = g_map.getGridInfoVoid();
659
660  while (it2.isNext())
661  {
662  auto key = it2.get();
663  auto key_g = g_map.getGKey(key);
664  key_g += pd.getKP1();
665
666  // To linearize must be positive
667  bool is_negative = false;
668  for (size_t i = 0 ; i < Sys_eqs::dims ; i++)
669  {
670  if (key_g.get(i) < 0)
671  is_negative = true;
672  }
673
674  if (is_negative == true)
675  {
676  ++it2;
677  continue;
678  }
679
680  // Carefull g map is extended, so the original (0,0) is shifted in g_map by
681
682  if (g_map.template get<0>(key) == 7)
683  {
684  int debug = 0;
685  debug++;
686  }
687
688  map_row[g_map.template get<0>(key)] = ginfo.LinId(key_g);
689
690  ++it2;
691  }
692
694
695  Vcluster & v_cl = *global_v_cluster;
696
697  openfpm::vector<triplet> & trpl = A.getMatrixTriplets();
698
699  auto it = it_d;
700  grid_sm<Sys_eqs::dims,void> gs = g_map.getGridInfoVoid();
701
702  std::unordered_map<long int,float> cols;
703
704  // resize b if needed
705  b.resize(Sys_eqs::nvar * g_map.size());
706
707  bool is_first = skip_first;
708
709  // iterate all the grid points
710  while (it.isNext())
711  {
712  if (is_first == true && v_cl.getProcessUnitID() == 0)
713  {
714  ++it;
715  is_first = false;
716  continue;
717  }
718  else
719  is_first = false;
720
721  // get the position
722  auto key = it.get();
723
724  // Calculate the non-zero colums
725  T::value(g_map,key,gs,spacing,cols,1.0);
726
728
729  auto g_calc_pos = g_map.getGKey(key);
730  g_calc_pos += pd.getKP1();
731
733
734  // create the triplet
735
736  for ( auto it = cols.begin(); it != cols.end(); ++it )
737  {
739  trpl.last().row() = g_map.template get<0>(key)*Sys_eqs::nvar + id;
740  trpl.last().col() = it->first;
741  trpl.last().value() = it->second;
742
744
745  auto ginfo = g_map.getGridInfoVoid();
746
747  size_t r = (trpl.last().row() / Sys_eqs::nvar);
748  size_t r_rest = (trpl.last().row() % Sys_eqs::nvar);
749  size_t c = (trpl.last().col() / Sys_eqs::nvar);
750  size_t c_rest = (trpl.last().col() % Sys_eqs::nvar);
751  double val = trpl.last().value();
752
753  // Transform
754
755  size_t rf = map_row[r] * 3 + r_rest;
756  size_t cf = map_row[c] * 3 + c_rest;
757
758  auto position_row = ginfo.InvLinId(rf / 3);
759  auto position_col = ginfo.InvLinId(cf / 3);
760
761  double valf = Al.getValue(rf,cf);
762
763  if (val != valf)
764  {
765  int debug = 0;
766  debug++;
767  }
768
770
771 // std::cout << "(" << trpl.last().row() << "," << trpl.last().col() << "," << trpl.last().value() << ")" << "\n";
772  }
773
774  b(g_map.template get<0>(key)*Sys_eqs::nvar + id) = num;
775
776  cols.clear();
777 // std::cout << "\n";
778
779  // if SE_CLASS1 is defined check the position
780 #ifdef SE_CLASS1
781 // T::position(key,gs,s_pos);
782 #endif
783
784  ++row;
785  ++row_b;
786  ++it;
787  }
788 }
789
790 typename Sys_eqs::SparseMatrix_type A;
791
797 typename Sys_eqs::SparseMatrix_type & getA()
798 {
799 #ifdef SE_CLASS1
800  consistency();
801 #endif
802  A.resize(g_map.size()*Sys_eqs::nvar,g_map.size()*Sys_eqs::nvar);
803
805
806 // A.save("debug_matrix_single_processor");
807
809
810  return A;
811
812 }
813
814
815 typename Sys_eqs::SparseMatrix_type A;
816
822 typename Sys_eqs::SparseMatrix_type & getA()
823 {
824 #ifdef SE_CLASS1
825  consistency();
826 #endif
827  A.resize(g_map.size()*Sys_eqs::nvar,g_map.size()*Sys_eqs::nvar);
828
830
831 // A.save("debug_matrix_single_processor");
832
834
835  return A;
836
837 }
838
839
845 typename Sys_eqs::Vector_type & getB()
846 {
847 #ifdef SE_CLASS1
848  consistency();
849 #endif
850
851  // size of the matrix
852 // B.resize(g_map.size()*Sys_eqs::nvar);
853
854  // copy the vector
855 // for (size_t i = 0; i < row_b; i++)
856 // B.insert(i,b.get(i));
857
858  return b;
859 }
860 };
861
862
863
864
870 {
871 /*
872 #ifdef SE_CLASS1
873
874  // No box must be unlinked
875  for (size_t i = 0 ; i < proc_int_box.size() ; i++)
876  {
877  for (size_t j = 0 ; j < proc_int_box.get(i).ibx.size() ; j++)
879
880  for (size_t j = 0 ; j < proc_int_box.get(i).ebx.size() ; j++)
882  }
883 #endif
884
885  // Get the number of near processors
886  for (size_t i = 0 ; i < proc_int_box.size() ; i++)
887  {
888  std::unordered_map<size_t,std::pair<size_t,size_t>> from_id_to_ibox;
889  std::unordered_map<size_t,std::pair<size_t,size_t>> from_id_to_ebox;
890
891  for (size_t j = 0 ; j < getProcessorNIGhost(i) ; j++)
892  {
893  std::pair<size_t,size_t> & ele = from_id_to_ibox[getProcessorIGhostId(i,j)];
894  ele.first = i;
895  ele.second = j;
896  }
897
898  for (size_t j = 0 ; j < getProcessorNEGhost(i) ; j++)
899  {
900  std::pair<size_t,size_t> & ele = from_id_to_ebox[getProcessorEGhostId(i,j)];
901
902  ele.first = i;
903  ele.second = j;
904  }
905
906  // iterate across all the ibox
907
908  for ( auto it = from_id_to_ibox.begin(); it != from_id_to_ibox.end(); ++it )
909  {
910  auto ite = from_id_to_ebox.find(it->first);
911
912  if(ite == from_id_to_ebox.end())
913  std::cerr << __FILE__ << ":" << __LINE__ << " error exist an internal ghost box that does not have an external ghost box" << std::endl;
914
915  if (ite->first != it->first)
916  std::cerr << __FILE__ << ":" << __LINE__ << " error exist an internal ghost box with inconsistent information about its origin" << std::endl;
917
920  }
921  }
922
923 #ifdef SE_CLASS1
924
925  // No box must be unlinked
926  for (size_t i = 0 ; i < proc_int_box.size() ; i++)
927  {
928  for (size_t j = 0 ; j < proc_int_box.get(i).ibx.size() ; j++)
929  {
931  std::cerr << __FILE__ << ":" << __LINE__ << " error detected unlinked internal ghost box" << std::endl;
932  }
933
934  for (size_t j = 0 ; j < proc_int_box.get(i).ebx.size() ; j++)
935  {
937  std::cerr << __FILE__ << ":" << __LINE__ << " error detected unlinked external ghost box" << std::endl;
938  }
939  }
940 #endif*/
941
942
943  /* for (size_t i = 0 ; i < this->getNNProcessors() ; i++)
944  {
945  for (size_t j = 0 ; j < this->getProcessorNIGhost(i) ; j++)
946  {
947  size_t id_i = this->getProcessorIGhostId(i,j);
949
951  return false;
952
954
955  if (id_i != id_e)
956  return false;
957  }
958  }*/
959 }
960
961
962
963
964
966
981 inline bool fix_box_ig(Box<dim,size_t> & bs, Box<dim,long int> & dom_i, const Box<dim,size_t> & bd, comb<dim> & cmb)
982 {
983  // Each dimension must match
984  for (size_t k = 0 ; k < dim ; k++)
985  {
986  size_t iw = bs.getHigh(k) - bs.getLow(k);
987  size_t ew = bd.getHigh(k) - bd.getLow(k);
988
989  if (iw != ew)
990  {
991  std::cout << "Fixing internal external" << std::endl;
992
993  Box<dim,size_t> & bst = bs;
994
995  if (cmb.c[k] == -1)
996  bst.setHigh(k,bd.getHigh(k) - (iw - ew));
997  else if (cmb.c[k] == 1)
998  bst.setLow(k,bs.getLow(k) + (iw - ew));
999  else
1000  return false;
1001
1002  // points in direction k of the domain
1003  long int dom_ext = dom_i.getHigh(k) - dom_i.getLow(k);
1004  // points in direction k of the internal ghost box
1005  long int ext_ibox = bst.getHigh(k) - bst.getLow(k);
1006
1007  // internal ghost box cannot be bigger than the domain
1008  // notify the failure in fixing
1009  if (dom_ext < ext_ibox)
1010  return false;
1011
1012  bs = bst;
1013  }
1014  }
1015
1016  return true;
1017 }
1018
1020
1021
1022 bool ret = fix_box_ig(bx_src,gdb_ext.get(i).Dbox,bx_dst,loc_eg_box.get(sub_id_dst).bid.get(k).cmb);
1023
1024 if (ret == false)
1025  std::cerr << "ERROR FAIL TO FIX " << std::endl;
1026
1027
1029
1030
1042 inline void fix_ie_g_box()
1043 {
1044  if (init_fix_ie_g_box == true) return;
1045
1046  comb<dim> zero;
1047  zero.zero();
1048
1049  // Here we collect all the external ghost box in the sector different from 0 that this processor has
1050
1052  openfpm::vector<size_t> prc_recv;
1053  openfpm::vector<size_t> sz_recv;
1054  openfpm::vector<openfpm::vector<Box_fix<dim>>> box_ext_send(dec.getNNProcessors());
1056
1057  // It contain the map g_id as key, and the pair, processor id, box-id
1058  std::unordered_map<long int,std::pair<long int,long int>> iglist;
1059
1060  // Here we create list of all the internal ghost box linked with an external ghost box
1061  // by periodicity
1062  for(size_t i = 0 ; i < dec.getNNProcessors() ; i++)
1063  {
1064  for (size_t j = 0 ; j < ig_box.get(i).bid.size() ; j++)
1065  {
1066  if (ig_box.get(i).bid.get(j).cmb != zero)
1067  {
1068  auto & ele = iglist[ig_box.get(i).bid.get(j).g_id];
1069  ele.first = i;
1070  ele.second = j;
1071  }
1072  }
1073  }
1074
1075  for(size_t i = 0 ; i < dec.getNNProcessors() ; i++)
1076  {
1077  for (size_t j = 0 ; j < eg_box.get(i).bid.size() ; j++)
1078  {
1079  if (eg_box.get(i).bid.get(j).cmb != zero)
1080  {
1082  box_ext_send.get(i).last().bx = eg_box.get(i).bid.get(j).l_e_box;
1083  box_ext_send.get(i).last().g_id = eg_box.get(i).bid.get(j).g_id;
1084  }
1085  }
1087  }
1088
1089  v_cl.SSendRecv(box_ext_send,box_ext_recv,prc,prc_recv,sz_recv);
1090
1091  // Received the external boxes we do fixation for each processor
1092  for (size_t i = 0 ; i < box_ext_recv.size() ; i++)
1093  {
1094  // For each received external ghost box
1095  for (size_t j = 0 ; j < box_ext_recv.get(i).size() ; j++)
1096  {
1098  size_t proc_id = dec.ProctoID(prc_recv.get(i));
1099
1100  auto it = g_id_to_internal_ghost_box.get(proc_id).find(box_ext_recv.get(i).get(j).g_id);
1101
1102 #ifdef SE_CLASS1
1103
1104  if (it == g_id_to_internal_ghost_box.get(proc_id).end())
1105  {
1106  std::cerr << __FILE__ << ":" << __LINE__ << " warning unlinked external ghost box" << std::endl;
1107  continue;
1108  }
1109
1110 #endif
1111
1113
1114  Box<dim,size_t> & box_i = ig_box.get(proc_id).bid.get(link).box;
1115
1116  // local Sub-domain from where this internal ghost box is calculated
1117  Box<dim,long int> & box_sub_i = gdb_ext.get(ig_box.get(proc_id).bid.get(link).sub).Dbox;
1118
1120
1121  // the fixing can fail
1122  // if it fail put the ig_box into a list
1123  // The fix can fail (for example) if the external ghost box require 7 point on x
1124  // but the domain has 6 point, in this case we cannot correct the internal ghost box
1125  bool ret = fix_box_ig(box_i,box_sub_i,box_ext_recv.get(i).get(j).bx,cmb);
1126
1127  if (ret == false)
1128  std::cerr << __FILE__ << ":" << __LINE__ << " and inconsistency between internal and external ghost boxes has been detected. The fix is not possible please change your ghost size (by a small amount) on the order of 10^-5 if you use float 10^-14 if you use double" << std::endl;
1129
1130  // Invalidate the ig_box in the list
1131  auto & ele = iglist[box_ext_recv.get(i).get(j).g_id];
1132  ele.first = -1;
1133  ele.second = -1;
1134  }
1135  }
1136
1137  // Here we check if all the internal ghost box has been explored
1138  // if one internal ghost box has not been explored, it been that, there is not
1139  // corresponding external ghost box on the other side. so we invalidate
1140
1141  for ( auto it = iglist.begin(); it != iglist.end(); ++it )
1142  {
1143  // If has not been explored invalidate, there is not external ghost
1144  if (it->second.first != -1)
1145  {
1146  size_t a = it->second.first;
1147  size_t b = it->second.second;
1148  ig_box.get(a).bid.get(b).box.invalidate();
1149  }
1150  }
1151 }
1152
1153
1155
1156 // Fix the exteenal and internal ghost boxes in ghost get
1157 fix_ie_g_box();
1158
1160
1161
1162 /* Point<dim,long int> p;
1163  p.get(0) = 0;
1164  p.get(1) = 81;
1165  p.get(2) = 79;
1166  if (ib.isInside(p))
1167  {
1168  int debug = 0;
1169  debug++;
1170  }
1171
1172  for (size_t i = 0 ; i < dim ; i++)
1173  {
1174  if (sub_domain.getLow(i) == ib_dom.getLow(i) &&
1175  (sub_domain_other.getHigh(i) == sub_domain.getLow(i) || cmb.c[i] == 1))
1176  {
1177  if (g.getHigh(i) != INVALID_GHOST && (ib.getHigh(i) - ib.getLow(i) + 1) > g.getHigh(i))
1178  {
1179  ib.setHigh(i,ib.getLow(i) + g.getHigh(i) - 1);
1180  }
1181  }
1182
1183  if (sub_domain.getHigh(i) == ib_dom.getHigh(i) &&
1184  (sub_domain_other.getLow(i) == sub_domain.getHigh(i) || cmb.c[i] == 1))
1185  {
1186  if (g.getLow(i) != -INVALID_GHOST && (ib.getHigh(i) - ib.getLow(i) + 1) > abs(g.getLow(i)))
1187  {
1188  ib.setLow(i, g.getHigh(i) - g.getLow(i) + 1);
1189  }
1190  }
1191  }
1192
1193  // This is a special case because a domain intersect itself by
1194  // periodicity
1195  if (sub_domain == sub_domain_other)
1196  {
1197  for (size_t i = 0 ; i < dim ; i++)
1198  {
1199  if (sub_domain.getLow(i) == ib_dom.getLow(i) &&
1200  sub_domain.getLow(i) == domain.getLow(i) &&
1201  sub_domain_other.getHigh(i) == domain.getHigh(i) &&
1202  cmb.c[i] == 1)
1203  {
1204  if (g.getHigh(i) != INVALID_GHOST && (ib.getHigh(i) - ib.getLow(i) + 1) > g.getHigh(i))
1205  {
1206  ib.setHigh(i,ib.getLow(i) + g.getHigh(i) - 1);
1207  }
1208  }
1209
1210  if (sub_domain.getHigh(i) == ib_dom.getHigh(i) &&
1211  sub_domain.getHigh(i) == domain.getHigh(i) &&
1212  sub_domain_other.getLow(i) == sub_domain.getHigh(i) &&
1213  cmb.c[i] == -1)
1214  {
1215  if (g.getLow(i) != -INVALID_GHOST && (ib.getHigh(i) - ib.getLow(i) + 1) > abs(g.getLow(i)))
1216  {
1217  ib.setLow(i, g.getHigh(i) - g.getLow(i) + 1);
1218  }
1219  }
1220  }
1221  }*/
1222
1224
1225 #endif /* GARGABE_HPP_ */
This class represent an N-dimensional box.
Definition: SpaceBox.hpp:26
static size_t getNumberOfElements_R(size_t d)
Get the number of Elements of dimension d.
Definition: HyperCube.hpp:67
size_t getProcessUnitID()
Get the process unit id.
static std::vector< comb< dim > > getCombinations_R(size_t d)
Definition: HyperCube.hpp:100
class that store Internal part external and border part of a dataset
grid_key_dx is the key to access any element in the grid
Definition: grid_key.hpp:18
Position of the element of dimension d in the hyper-cube of dimension dim.
Definition: comb.hpp:34
This represent a sub-hyper-cube of an hyper-cube like a face or an edge of a cube.
Definition: HyperCube.hpp:8
__device__ __host__ T getLow(int i) const
get the i-coordinate of the low bound interval of the box
Definition: Box.hpp:556
boost::shared_ptr< T > ext
__device__ __host__ index_type get(index_type i) const
Get the i index.
Definition: grid_key.hpp:503
This class implement the point shape in an N-dimensional space.
Definition: Point.hpp:27
size_t size()
Stub size.
Definition: map_vector.hpp:211
T getVolumeKey() const
Get the volume spanned by the Box P1 and P2 interpreted as grid key.
Definition: Box.hpp:1351
Implementation of VCluster class.
Definition: VCluster.hpp:58
__device__ __host__ void setHigh(int i, T val)
set the high interval of the box
Definition: Box.hpp:544
void execute()
Execute all the requests.
This class decompose a space into sub-sub-domains and distribute them across processors.
Sparse Matrix implementation.
This is a distributed grid.
__device__ __host__ void setLow(int i, T val)
set the low interval of the box
Definition: Box.hpp:533
boost::shared_ptr< T > inte
Distributed grid iterator.
boost::shared_ptr< T > bord
Border part of the data.
void zero()
Set all the elements to zero.
Definition: comb.hpp:83
This class represent an N-dimensional box.
Definition: Box.hpp:60
This class calculate elements of the hyper-cube.
Definition: HyperCube.hpp:57
Implementation of 1-D std::vector like structure.
Definition: map_vector.hpp:202
__device__ __host__ T getHigh(int i) const
get the high interval of the box
Definition: Box.hpp:567
bool SSendRecv(openfpm::vector< T > &send, S &recv, openfpm::vector< size_t > &prc_send, openfpm::vector< size_t > &prc_recv, openfpm::vector< size_t > &sz_recv, size_t opt=NONE)
Semantic Send and receive, send the data to processors and receive from the other processors.
Definition: VCluster.hpp:797
char c[dim]
Array that store the combination.
Definition: comb.hpp:37