Allink  v0.1
Cubo.mirror.cpp
1 #include <Cubo.h>
2 #define MASK_X 0
3 #define MASK_Y 4
4 #define MASK_Z 8
5 #define GET_MASK_X 0x00f
6 #define GET_MASK_Y 0x0f0
7 #define GET_MASK_Z 0xf00
8 #define DOM_CENTER 0
9 #define DOM_LEFT 1
10 #define DOM_RIGHT 2
11 #define DOM_SET_POS(n,p,m) ( (n)|=((p)<<(m)) )
12 #define DOM_GET_POS(n,m,g) ( ((n)&(g))>>(m) )
13 //FIXME
14 #define DOM_IF_POS(n,m,d) ( ((n)>>(m))==(d) )
15 //-----------------BASICS-----------------------------
16 void DomDecBasics::SigErr(int IsWrong,const char * s, ...){
17 #ifdef DEBUG
18  if(IsWrong){
19  va_list args;
20  va_start(args, s);
21  fprintf(stderr, "Error: ");
22  vfprintf(stderr, s, args);
23  fprintf(stderr, "exiting the program\n");
24  va_end(args);
25  abort();
26  }
27 #else
28  return;
29 #endif
30 }
31 DomDecBasics::DomDecBasics(double EdgeExt[3],int NPartExt,double CutOffExt){
32  IfMirror = 1;
33  IfBorder = 0;
34  double CellSize = IfBorder ? CutOffExt : 2.*CutOffExt;
35  for(int d=0;d<3;d++){
36  Edge[d] = EdgeExt[d];
37  InvEdge[d] = 1./Edge[d];
38  NSect[d] = (int)(Edge[d]/CellSize) + 2*IfMirror;
39  SigErr(CutOffExt > 2.*Edge[d],"Insufficient number of cells \n");
40  }
41  NPart = 0;
42  SetCutOff(CutOffExt);
43  Mod10[0] = 1;
44  Mod10[1] = NSect[0];
45  Mod10[2] = NSect[1]*Mod10[1];
46  NCell = NSect[0]*NSect[1]*NSect[2];
47 }
48 int DomDecBasics::SetCutOff(double CutOffExt){
49  for(int d=0;d<3;d++){
50  if(CutOffExt > Edge[d]/(double)NSect[d]){
51  printf("CutOff larger than the cell size\n");
52  return 1;
53  }
54  }
55  CutOff = CutOffExt;
56  return 0;
57 }
59  return NPart;
60 }
62  return NCell;
63 }
64 int DomDecBasics::GetCell(double *Pos,int *NeiList){
65  int c = pCella(Pos);
66  return GetCellCh(c,NeiList);
67 }
68 int DomDecBasics::pCella(const double Pos[3]){
69  int ris[3];
70  for(int d=0;d<3;d++){
71  ris[d] = (int)(Pos[d]*InvEdge[d]*NSect[d]);
72  }
73  return pCella(ris);
74 }
75 int DomDecBasics::pCella(const double *Pos,int *cd){
76  for(int d=0;d<3;d++){
77  cd[d] = (int)(Pos[d]*InvEdge[d]*NSect[d]);
78  }
79  return pCella(cd);
80 }
81 int DomDecBasics::pCella(const int c[3]){
82  //+1 for the mirror notation
83  int risp = (c[0]+IfMirror)*Mod10[0]+(c[1]+IfMirror)*Mod10[1]+(c[2]+IfMirror)*Mod10[2];
84  SigErr(risp >= NCell,"Probably the particle position is not a number or not in the box\n");
85  return risp;
86 }
87 int DomDecBasics::GetCellChMirror(int c[3],int *NeiList){
88  int NeiMod[3];
89  int NeiCell[3];
90  int NNei = 0;
91  for(int cx=-1;cx<=1;cx++){
92  NeiCell[0] = NeiMod[0] + cx;
93  for(int cy=-1;cy<=1;cy++){
94  NeiCell[1] = NeiMod[1] + cy;
95  for(int cz=-1;cz<=1;cz++){
96  NeiList[NNei++] = NeiCell[0]*Mod10[0] + NeiCell[1]*Mod10[1] + NeiCell[2]*Mod10[2];
97  }
98  }
99  }
100  return NNei;
101 }
102 int DomDecBasics::GetCellCh(int c,int *NeiList){
103  int cd[3];
104  cd[2] = (int)(c/(double)Mod10[2]);
105  cd[1] = (int)((c-cd[2]*Mod10[2])/(double)Mod10[1]);
106  cd[0] = (int)(c-cd[1]*Mod10[1]-cd[2]*Mod10[2]);
107  GetCellCh(cd,NeiList);
108 }
109 int DomDecBasics::GetCellCh(int *cd,int *NeiList){
110  int NeiMod[3] = {cd[0],cd[1],cd[2]};
111  int NeiCell[3];
112  int NNei = 0;
113  for(int cx=-1;cx<=1;cx++){
114  NeiCell[0] = NeiMod[0] + cx;
115  if(NeiCell[0] >= NSect[0]) NeiCell[0] -= NSect[0];
116  if(NeiCell[0] < 0) NeiCell[0] += NSect[0];
117  for(int cy=-1;cy<=1;cy++){
118  NeiCell[1] = NeiMod[1] + cy;
119  if(NeiCell[1] >= NSect[1]) NeiCell[1] -= NSect[1];
120  if(NeiCell[1] < 0) NeiCell[1] += NSect[1];
121  for(int cz=-1;cz<=1;cz++){
122  NeiCell[2] = NeiMod[2] + cz;
123  if(NeiCell[2] >= NSect[2]) NeiCell[2] -= NSect[2];
124  if(NeiCell[2] < 0) NeiCell[2] += NSect[2];
125  NeiList[NNei++] = NeiCell[0]*Mod10[0] + NeiCell[1]*Mod10[1] + NeiCell[2]*Mod10[2];
126  }
127  }
128  }
129  return NNei;
130 }
131 int DomDecBasics::GetCellCoord(double *Pos,int *NeiList){
132  int c = pCella(Pos);
133  int Coord = GetCoorNumb(Pos);
134  return GetCellCoord(c,Coord,NeiList);
135 }
136 int DomDecBasics::GetCellCoord(int c,int Coord,int *NeiList){
137  SigErr(c < 0,"Probably the particle position is not a number\n");
138  int NNei = 0;
139  int nNei[3];
140  int NeiCell[3];
141  int Mask[3] = {MASK_X,MASK_Y,MASK_Z};
142  int GetMask[3] = {GET_MASK_X,GET_MASK_Y,GET_MASK_Z};
143  int Move[3] = {0,0,0};
144  NeiCell[2] = (int)(c/(double)Mod10[2]);
145  NeiCell[1] = (int)((c-NeiCell[2]*Mod10[2])/(double)Mod10[1]);
146  NeiCell[0] = (int)(c-NeiCell[1]*Mod10[1]-NeiCell[2]*Mod10[2]);
147  for(int d=0;d<3;d++){
148  int d1 = (d+1)%3;
149  int d2 = (d+2)%3;
150  if( DOM_GET_POS(Coord,Mask[d],GetMask[d]) == DOM_LEFT ){
151  Move[d] = DOM_LEFT;
152  int nc = NeiCell[d] - 1;
153  if(nc < 0) nc = NSect[d]-1;
154  NeiList[NNei++] = nc*Mod10[d] + NeiCell[d1]*Mod10[d1] + NeiCell[d2]*Mod10[d2];
155  }
156  else if( DOM_GET_POS(Coord,Mask[d],GetMask[d]) == DOM_RIGHT ){
157  Move[d] = DOM_RIGHT;
158  int nc = NeiCell[d] + 1;
159  if(nc >= NSect[d]) nc = 0;
160  NeiList[NNei++] = nc*Mod10[d] + NeiCell[d1]*Mod10[d1] + NeiCell[d2]*Mod10[d2];
161  }
162  }
163  if(NNei > 1){
164  for(int d=0;d<3;d++){
165  nNei[d] = NeiCell[d];
166  if(Move[d] == DOM_LEFT){
167  nNei[d] = NeiCell[d] - 1;
168  if(nNei[d] < 0) nNei[d] = NSect[d] - 1;
169  }
170  if(Move[d] == DOM_RIGHT){
171  nNei[d] = NeiCell[d] + 1;
172  if(nNei[d] >= NSect[d]) nNei[d] = 0;
173  }
174  }
175  NeiList[NNei++] = nNei[0]*Mod10[0] + nNei[1]*Mod10[1] + nNei[2]*Mod10[2];
176  }
177  if(NNei == 4){
178  for(int dd=0;dd<3;dd++){
179  for(int d=0;d<3;d++){
180  nNei[d] = NeiCell[d];
181  if(d == dd) continue;
182  if(Move[d] == DOM_LEFT){
183  nNei[d] = NeiCell[d] - 1;
184  if(nNei[d] < 0) nNei[d] = NSect[d] - 1;
185  }
186  if(Move[d] == DOM_RIGHT){
187  nNei[d] = NeiCell[d] + 1;
188  if(nNei[d] >= NSect[d]) nNei[d] = 0;
189  }
190  }
191  NeiList[NNei++] = nNei[0]*Mod10[0] + nNei[1]*Mod10[1] + nNei[2]*Mod10[2];
192  }
193  }
194  NeiList[NNei++] = c;
195  // for(int n=0;n<NNei;n++) assert(NeiList[n] < NCell);
196  // for(int n=0;n<NNei;n++) assert(NeiList[n] >= 0);
197  return NNei;
198 }
199 int DomDecBasics::GetCoorNumb(double *Pos){
200  double CellLen[3];
201  int Border[3];
202  double BorderR[3];
203  int Mask[3] = {MASK_X,MASK_Y,MASK_Z};
204  int Coord = 0;
205  for(int d=0;d<3;d++){
206  Border[d] = (int)(Pos[d]/Edge[d]*NSect[d]);
207  CellLen[d] = (Edge[d]/(double)NSect[d]);
208  if(Pos[d] - Border[d]*CellLen[d] < CutOff){
209  DOM_SET_POS(Coord,DOM_LEFT,Mask[d]);
210  }
211  if(-Pos[d] + (Border[d]+1)*CellLen[d] < CutOff){
212  DOM_SET_POS(Coord,DOM_RIGHT,Mask[d]);
213  }
214  }
215  return Coord;
216 }
217 //-------------------LINKED-LIST-------------------------------
218 void DdLinkedList::Erase(){
219  for(int c=0;c<NCell;c++){
220  Cella[c].NPart = 0;
221  Cella[c].Last = -1;
222  }
223  // for(int p=0;p<NPart;p++){
224  // Pc[p].Cell = 0;
225  // Pc[p].Next = -1;
226  // Pc[p].Prev = -2;
227  // }
228  NPart = 0;
229 }
230 int DdLinkedList::SetCoorNumb(double *Pos,int p){
231  Pc[p].Coord = GetCoorNumb(Pos);
232  return Pc[p].Coord;
233 }
234 void DdLinkedList::SetCounters(int c1){
235  for(int c=0;c<NCell;c++){
236  Cella[c].Curr1 = Cella[c].First;
237  Cella[c].Curr2 = Cella[c].First;
238  if(Cella[c].NPart > 1)
239  Cella[c].Curr2 = Pc[Cella[c].Curr1].Next;
240  }
241 }
242 void DdLinkedList::AddPart(const int p,double *Pos){
243  int cd[3];
244  pCella(Pos,cd);
245  SetCoorNumb(Pos,p);
246  int c = pCella(cd);
247  int pl = Cella[c].Last;
248  SigErr(pl == p,"The particle %d is already the last in the cell %d\n",p,c);
249  SigErr(p >= NAllocP,"More particle than allocated\n");
250  SigErr(Cella[c].NPart >= NAllocP,"More particle than allocated\n");
251  Cella[c].NPart++;
252  NPart++;
253  SigErr(pl >= NAllocP,"More particle than allocated\n");
254  Cella[c].Last = p;
255  for(int d=0;d<3;d++){
256  Pc[p].Cell[d] = cd[d];
257  Pc[p].Pos[d] = Pos[d];
258  }
259  AddPartMirror(p);
260  Pc[p].Next = -1;
261  if(Cella[c].NPart == 1){
262  Cella[c].First = p;
263  Pc[p].Prev = -2;
264  return ;
265  }
266  if(pl >= 0) Pc[pl].Next = p;
267  Pc[p].Prev = pl;
268  //printf("%d - %d - %d c) %d last %d coord %x cell %d\n",pp,p,pn,c,Cella[c].Last,Pc[p].Coord,Pc[p].Cell);
269 }
270 void DdLinkedList::CheckMirror(int p){
271  int NCellMirror = 0;
272  for(int d=0;d<3;d++){
273  double Pos[3] = {Pc[p].Pos[0],Pc[p].Pos[1],Pc[p].Pos[2]};
274  int cd[3] = {Pc[p].Cell[0],Pc[p].Cell[1],Pc[p].Cell[2]};
275  if(Pc[p].Pos[d] - CutOff < 0.){
276  Pos[d] = Pc[p].Pos[d] - CutOff;
277  cd[d] -= 1;
278  AddPart(p,cd,Pos);
279  }
280  else if(Edge[d] - Pc[p].Pos[d] < CutOff){
281  Pos[d] = Pc[p].Pos[d] + CutOff;
282  cd[d] += 1;
283  AddPart(p,cd,Pos);
284  }
285  }
286 }
287 void DdLinkedList::AddMirrorPart(const int p,int *cd,double *Pos){
288  SetCoorNumb(Pos,p);
289  int c = pCella(cd);
290  int pl = Cella[c].Last;
291  SigErr(pl == p,"The particle %d is already the last in the cell %d\n",p,c);
292  SigErr(p >= NAllocP,"More particle than allocated\n");
293  SigErr(Cella[c].NPart >= NAllocP,"More particle than allocated\n");
294  Cella[c].NPart++;
295  NPart++;
296  SigErr(pl >= NAllocP,"More particle than allocated\n");
297  Cella[c].Last = p;
298  for(int d=0;d<3;d++){
299  Pc[p].Cell[d] = cd[d];
300  Pc[p].Pos[d] = Pos[d];
301  }
302  Pc[p].Next = -1;
303  if(Cella[c].NPart == 1){
304  Cella[c].First = p;
305  Pc[p].Prev = -2;
306  return ;
307  }
308  if(pl >= 0) Pc[pl].Next = p;
309  Pc[p].Prev = pl;
310 }
311 void DdLinkedList::RemPart(const int p,const int c){
312  SigErr(p >= NAllocP,"More particles than allocated\n");
313  //assert(Cella[c].NPart > 0);
314  SigErr(Cella[c].NPart <= 0,"Cannot remove, the cell is already empty\n");
315  SigErr(c < 0,"The cell %d does not exist\n",c);
316  Cella[c].NPart--;
317  NPart--;
318  if(Cella[c].NPart == 0){
319  Cella[c].First = -2;
320  Cella[c].Last = -1;
321  return ;
322  }
323  int pp = Pc[p].Prev;
324  int pn = Pc[p].Next;
325  if(pp == -2) Cella[c].First = pn;
326  else Pc[pp].Next = pn;
327  if(pn == -1) Cella[c].Last = pp;
328  else Pc[pn].Prev = pp;
329  for(int d=0;d<3;d++)
330  Pc[p].Cell[d] = -1;
331  //printf("%d - %d - %d c) %d last %d coord %x cell %d\n",pp,p,pn,c,Cella[c].Last,Pc[p].Coord,Pc[p].Cell);
332  // assert(pp != -2 && Pc[pp].Next != pp);
333 }
334 void DdLinkedList::SwapPart(int p1,double *Pos1,int p2,double *Pos2){
335  if(p1==p2) return ;
336  DOMAIN_PART Pn;
337  int c1 = pCella(Pc[p1].Cell);
338  int c2 = pCella(Pc[p2].Cell);
339  if(c1!=c2){
340  if(Cella[c2].First == p2)
341  Cella[c2].First = p1;
342  if(Cella[c2].Last == p2)
343  Cella[c2].Last = p1;
344  if(Cella[c1].First == p1)
345  Cella[c1].First = p2;
346  if(Cella[c1].Last == p1)
347  Cella[c1].Last = p2;
348  }
349  if(c1==c2){
350  if(Cella[c1].First == p2)
351  Cella[c1].First = p1;
352  else if(Cella[c1].First == p1)
353  Cella[c1].First = p2;
354  if(Cella[c1].Last == p2)
355  Cella[c1].Last = p1;
356  else if(Cella[c1].Last == p1)
357  Cella[c1].Last = p2;
358  }
359  int p1p = Pc[p1].Prev;
360  int p2p = Pc[p2].Prev;
361  int p1n = Pc[p1].Next;
362  int p2n = Pc[p2].Next;
363  if(p1p >= 0) Pc[p1p].Next = p2;
364  if(p1n >= 0) Pc[p1n].Prev = p2;
365  if(p2p >= 0) Pc[p2p].Next = p1;
366  if(p2n >= 0) Pc[p2n].Prev = p1;
367  Pn.Prev = Pc[p1].Prev;
368  Pn.Next = Pc[p1].Next;
369  Pn.Coord = Pc[p1].Coord;
370  for(int d=0;d<3;d++){
371  Pn.Cell[d] = Pc[p1].Cell[d];
372  Pc[p1].Cell[d] = Pc[p2].Cell[d];
373  Pc[p2].Cell[d] = Pn.Cell[d];
374  }
375  Pc[p1].Prev = Pc[p2].Prev;
376  Pc[p1].Next = Pc[p2].Next;
377  Pc[p1].Coord = Pc[p2].Coord;
378  Pc[p2].Prev = Pn.Prev;
379  Pc[p2].Next = Pn.Next;
380  Pc[p2].Coord = Pn.Coord;
381 }
382 void DdLinkedList::RemPart(const int p,double *Pos){
383  int c = pCella(Pos);
384  RemPart(p,c);
385 }
386 void DdLinkedList::RemPart(const int p){
387  int c = pCella(Pc[p].Cell);
388  SigErr(c < 0,"The cell %d does not exist\n",c);
389  return RemPart(p,c);
390 }
391 void DdLinkedList::MovePart(const int p,double *NewPos){
392  int cd[3];
393  int cn = pCella(NewPos,cd);
394  int co = pCell(p);
395  SetCoorNumb(NewPos,p);
396  if(cn == co) return;
397  RemPart(p,co);
398  AddPart(p,cd);
399 }
400 void DdLinkedList::MovePart(const int p,double *OldPos,double *NewPos){
401  int cd[3];
402  int cn = pCella(NewPos,cd);
403  int co = pCella(OldPos);
404  SetCoorNumb(NewPos,p);
405  if(cn == co) return ;
406  RemPart(p,co);
407  AddPart(p,cd);
408 }
409 int DdLinkedList::ItCell(const int c){
410  SigErr(c >= NCell,"The cell %d does not exist\n",c);
411  SigErr(Cella[c].Curr1 >= NAllocP,"Poiting to a particle %d over the number of allocated particles\n",Cella[c].Curr1,NAllocP);
412  return Cella[c].Curr1;
413 }
414 int DdLinkedList::IfItCell(const int c){
415  if(Cella[c].Curr1 < 0 || Cella[c].NPart == 0){
416  Cella[c].Curr1 = Cella[c].First;
417  return 0;
418  }
419  return 1;
420 }
421 void DdLinkedList::IncrCurr(const int c){
422  Cella[c].Curr1 = Pc[Cella[c].Curr1].Next;
423 }
424 void DdLinkedList::PrintCell(const int c){
425  for(SetCounters(c);IfItCell(c);IncrCurr(c)){
426  int p = ItCell(c);
427  printf("%d) # %d %d_%d_%d %x\n",c,Cella[c].NPart,Pc[p].Prev,ItCell(c),Pc[p].Next,Pc[p].Coord);
428  }
429 }
431  for(int c=0;c<NCell;c++)
432  PrintCell(c);
433 }
435  //It++;
436  return *this;
437 }
439 
440 
441 }
442 void DdLinkedList::Couple(const int c,int *p1,int *p2){
443  SigErr(Cella[c].Curr1 >= NAllocP,"Poiting to a particle %d over the number of allocated particles\n",Cella[c].Curr1,NAllocP);
444  SigErr(Cella[c].Curr2 >= NAllocP,"Poiting to a particle %d over the number of allocated particles\n",Cella[c].Curr2,NAllocP);
445  *p1 = Cella[c].Curr1;
446  *p2 = Cella[c].Curr2;
447 }
448 int DdLinkedList::IfItCouple(const int c){
449  if(Cella[c].NPart < 2) return 0;
450  if(Cella[c].Curr1 == -1){
451  Cella[c].Curr1 = Cella[c].First;
452  Cella[c].Curr2 = Pc[Cella[c].First].Next;
453  return 0;
454  }
455  return 1;
456 }
457 void DdLinkedList::IncrCurrList(const int c){
458  int p2 = Cella[c].Curr2;
459  Cella[c].Curr2 = Pc[p2].Next;
460  if(Pc[p2].Next == -1){
461  int p1 = Cella[c].Curr1;
462  Cella[c].Curr1 = Pc[p1].Next;
463  Cella[c].Curr2 = Pc[Pc[p1].Next].Next;
464  if(Cella[c].Curr2 == -1){
465  Cella[c].Curr1 = -1;
466  }
467  }
468  SigErr(Cella[c].Curr2 >= NAllocP,"Poiting to a particle %d over the number of allocated particles\n",Cella[c].Curr2,NAllocP);
469 }
471  return *this;
472 }
473 void operator+(DomCell &Dc){
474 }
475 DomPart::DomPart(int ExtNPart){
476  NPart = ExtNPart;
477  Next = new int[NPart];
478  Prev = new int[NPart];
479 }
480 int DomPart::operator[](int col){
481  return Next[col];
482 }
483 //-------------------------------DdArray------------------------
487 void DdArray::Erase(){
488  for(int c=0;c<NCell;c++){
489  Cella[c].Part.clear();
490  }
491 }
492 void DdArray::AddPart(int p,double *Pos){
493  int c = pCella(Pos);
494  Cella[c].Part.push_back(p);
495  NPart++;
496 }
497 void DdArray::RemPart(int p,double *Pos){
498  int c = pCella(Pos);
499  Cella[c].Part.remove(p);
500  NPart--;
501 }
502 void DdArray::MovePart(int p,double *OldPos,double *NewPos){
503  int c1 = pCella(OldPos);
504  int c2 = pCella(NewPos);
505  if(c1 == c2) return;
506  Cella[c1].Part.remove(p);
507  Cella[c2].Part.push_back(p);
508 }
509 void DdArray::SwapPart(int p1,double *Pos1,int p2,double *Pos2){
510  int c1 = pCella(Pos1);
511  int c2 = pCella(Pos2);
512  if(c1 == c2) return;
513  Cella[c1].Part.remove(p1);
514  Cella[c1].Part.push_back(p2);
515  Cella[c2].Part.remove(p2);
516  Cella[c2].Part.push_back(p2);
517 }
518 void DdArray::SetCounters(int c){
519  //if(Cella[c].Part.size() == 0) return;
520  NCurr = Cella[c].Part.begin();
521  NCurr2 = Cella[c].Part.begin();
522  if(Cella[c].Part.size() > 1)
523  ++NCurr2;
524 };
525 int DdArray::IfItCell(int c){
526  if(NCurr == Cella[c].Part.end()) return 0;
527  return 1;
528 };
529 void DdArray::IncrCurr(int c){
530  ++NCurr;
531 }
532 int DdArray::ItCell(int c){
533  return *NCurr;
534 }
535 void DdArray::Couple(const int c,int *p1,int *p2){
536  *p1 = *NCurr;
537  *p2 = *NCurr2;
538 }
539 int DdArray::IfItCouple(const int c){
540  if(NCurr == Cella[c].Part.end())
541  return 0;
542  return 1;
543 }
544 void DdArray::IncrCurrList(const int c){
545  ++NCurr2;
546  if(NCurr2 == Cella[c].Part.end()){
547  NCurr2 = NCurr;
548  ++NCurr2;
549  ++NCurr;
550  }
551 }
552 void DdArray::PrintCell(const int c){
553  for(SetCounters(c);IfItCell(c);IncrCurr(c)){
554  int p = ItCell(c);
555  printf("%d) # %d %d\n",c,Cella[c].Part.size(),ItCell(c));
556  }
557 }
558 void DdArray::PrintCells(){
559  for(int c=0;c<NCell;c++)
560  PrintCell(c);
561 }
562 
563 
564 //--------------------------------DdFixedSize----------------------------
568 void DdFixedSize::Erase(){
569  for(int c=0;c<NCell;c++){
570  for(int p=0;p<Cella[c].NPart;p++){
571  Cella[c].Part[p] = -1;
572  }
573  Cella[c].NPart = 0;
574  }
575 }
576 void DdFixedSize::AddPart(int p,double *Pos){
577  int c = pCella(Pos);
578  AddPart(p,c);
579 }
580 void DdFixedSize::AddPart(int p,int c){
581  //printf("Adding %d in %d\n",p,c);
582  if(Cella[c].NPart > NPCell){
583  printf("Maximum allocated size in the domain decomposition (%d) reached, terminating\n",NPCell);
584  assert(Cella[c].NPart < NPCell);
585  }
586  Cella[c].Part[Cella[c].NPart++] = p;
587  NPart++;
588 }
589 void DdFixedSize::RemPart(int p1,double *Pos){
590  int c = pCella(Pos);
591  RemPart(p1,c);
592 }
593 void DdFixedSize::RemPart(int p1,int c){
594  //PrintCell(c);
595  //printf("removing %d from %d\n",p1,c);
596  int FoundPartInCell = 0;
597  int Temp = 0;
598  for(int p=0;p<Cella[c].NPart;p++){
599  if(p1==Cella[c].Part[p]){
600  for(int pp=p;pp<Cella[c].NPart-1;pp++){
601  Temp = Cella[c].Part[pp];
602  Cella[c].Part[pp] = Cella[c].Part[pp+1];
603  Cella[c].Part[pp+1] = Temp;
604  }
605  Cella[c].NPart--;
606  NPart--;
607  FoundPartInCell = 1;
608  break;
609  }
610  }
611  //PrintCell(c);
612  assert(FoundPartInCell);
613 }
614 void DdFixedSize::MovePart(int p,double *OldPos,double *NewPos){
615  int c1 = pCella(OldPos);
616  int c2 = pCella(NewPos);
617  if(c1 == c2) return;
618  RemPart(p,c1);
619  AddPart(p,c2);
620 }
621 void DdFixedSize::SwapPart(int p1,double *Pos1,int p2,double *Pos2){
622  int c1 = pCella(Pos1);
623  int c2 = pCella(Pos2);
624  if(c1 == c2) return;
625  RemPart(p1,c1);
626  RemPart(p2,c2);
627  AddPart(p1,c2);
628  AddPart(p2,c1);
629 }
630 void DdFixedSize::SetCounters(int c){
631  //if(Cella[c].Part.size() == 0) return;
632  NCurr = 0;
633  NCurr2 = 0;
634  if(Cella[c].NPart > 0) NCurr2++;
635 };
636 int DdFixedSize::IfItCell(int c){
637  if(NCurr == Cella[c].NPart) return 0;
638  return 1;
639 };
640 void DdFixedSize::IncrCurr(int c){
641  NCurr++;
642 }
643 int DdFixedSize::ItCell(int c){
644  return Cella[c].Part[NCurr];
645 }
646 void DdFixedSize::Couple(const int c,int *p1,int *p2){
647  // printf("%d %d %d\n",NCurr,NCurr2,Cella[c].NPart);
648  *p1 = Cella[c].Part[NCurr];
649  *p2 = Cella[c].Part[NCurr2];
650 }
651 int DdFixedSize::IfItCouple(const int c){
652  if(NCurr >= Cella[c].NPart)
653  return 0;
654  return 1;
655 }
656 void DdFixedSize::IncrCurrList(const int c){
657  NCurr2++;
658  if(NCurr2 >= Cella[c].NPart){
659  NCurr2 = NCurr;
660  NCurr2++;
661  NCurr++;
662  }
663 }
664 void DdFixedSize::PrintCell(const int c){
665  for(SetCounters(c);IfItCell(c);IncrCurr(c)){
666  int p = ItCell(c);
667  printf("%d) # %d %d\n",c,Cella[c].NPart,ItCell(c));
668  }
669 }
671  for(int c=0;c<NCell;c++)
672  PrintCell(c);
673 }
674 //-----------------------------NeiVertex------------------------
678 NeiVertex::NeiVertex(int NGridExt){
679  NGrid = 2*NGridExt-1;
680  int NVertex = NGrid*NGrid*NGrid;
681  Vertex = new VERTEX[NVertex];
682  vPos = new PPOS[NVertex];
683 }
685  delete [] Vertex;
686  delete [] vPos;
687 }
688 int NeiVertex::GetVertex(double *Pos){
689  int s[3];
690  for(int d=0;d<3;d++){
691  s[d] = (int)(Pos[d]*NGrid);
692  }
693  return (s[0]*NGrid+s[1])*NGrid+s[2];
694 }
695 void NeiVertex::Add(double *Pos,int t){
696  int v = GetVertex(Pos);
697  Vertex[v].v.push_back(t);
698  for(int d=0;d<3;d++)
699  vPos[v].Pos[d] = Pos[d];
700 }
701 void NeiVertex::Add(int v,int t,double *Pos){
702  Vertex[v].v.push_back(t);
703  for(int d=0;d<3;d++)
704  vPos[v].Pos[d] = Pos[d];
705 }
706 void NeiVertex::Rem(int v,int t){
707  Vertex[v].v.remove(t);
708 }
709 void NeiVertex::SetCounters(int v){
710  tCurr = Vertex[v].v.begin();
711 }
712 int NeiVertex::IfItCell(int v){
713  if(tCurr == Vertex[v].v.end()) return 0;
714  return 1;
715 };
716 void NeiVertex::IncrCurr(int v){
717  ++tCurr;
718 }
719 int NeiVertex::VertCurr(int v){
720  return *tCurr;
721 }
722 void NeiVertex::PosVertex(int v,double *Pos){
723  for(int d=0;d<3;d++)
724  Pos[d] = vPos[v].Pos[d];
725 }
726 
727 
728 //-----------------------------
int Coord
Coordination number.
Definition: Cubo.h:134
void AddPart(const int p, double *Pos)
Add a particle to the cell c.
Definition: Cubo.cpp:276
int IfItCouple(const int c)
Tell when the loop is over.
Definition: Cubo.cpp:701
void SigErr(int Condition, const char *s,...)
Handle errors.
Definition: Cubo.cpp:24
void PrintCells()
Print the the particles in the cells.
Definition: Cubo.cpp:446
void MovePart(int p, double *OldPos, double *NewPos)
Move a part.
Definition: Cubo.cpp:793
void Couple(const int c, int *p1, int *p2)
Returns the pointed couple of particle.
Definition: Cubo.cpp:825
int NCell
Number of cells.
Definition: Cubo.h:57
int VertCurr(int v)
Current vertex iterator for the vertex v.
Definition: Cubo.cpp:1109
void RemPart(int p, double *Pos)
Remove a part from the cell.
Definition: Cubo.cpp:768
int GetCellCoord(int c, int Coord, int *NeiList)
Get the number of the neighbouring cells for the particles close to a wall.
Definition: Cubo.cpp:121
int SetCutOff(double CutOffExt)
Set the CutOff.
Definition: Cubo.cpp:40
void Erase()
Erase the content of the cells.
Definition: Cubo.cpp:747
void SetCounters(int c)
Set the counters to the initial position.
Definition: Cubo.cpp:413
int IfItCell(int c)
End the iteration in the cell.
Definition: Cubo.cpp:815
void Erase()
Erase the pairlist.
Definition: Cubo.cpp:258
void IncrCurrList(const int c)
Increment the pointers.
Definition: Cubo.cpp:706
void PrintCell(const int c)
Print the content of one cell.
Definition: Cubo.cpp:843
int NSect[3]
Number of cells per direction.
Definition: Cubo.h:59
int IfItCouple(const int c)
Stop the loop and set the counter to zero.
Definition: Cubo.cpp:572
DomPart & operator++()
Increase the iterator, deactivated.
Definition: Cubo.cpp:616
int GetCellCh(int c, int *NeiList)
Get the number of the neighbouring cells.
Definition: Cubo.cpp:64
void PrintCell(const int c)
Print the particles in a cell.
Definition: Cubo.cpp:439
~NeiVertex()
Delete the class.
Definition: Cubo.cpp:1022
void AddPart(int p, double *Pos)
Add a part to the cell.
Definition: Cubo.cpp:654
void SetCounters(int c)
Set the counters to the initial pointer.
Definition: Cubo.cpp:809
NeiVertex(int NTriaExt, int NvPtExt, int NGridExt, double *EdgeExt)
Alloc NVertices.
Definition: Cubo.cpp:1007
int IfItCell(int v)
End of the counters.
Definition: Cubo.cpp:1102
void RemPart(const int p, double *Pos)
Remove a particle form the cell c.
Definition: Cubo.cpp:307
int NPart
Total number of particle.
Definition: Cubo.h:71
void SwapPart(int p1, double *Pos1, int p2, double *Pos2)
Swap two parts.
Definition: Cubo.cpp:800
int ItCell(int c)
Value of the current particle in the cell c.
Definition: Cubo.cpp:694
DomCell & operator=(const DomCell &Dc)
Return the cell class.
Definition: Cubo.cpp:619
void PrintCells()
Print the content of every cell.
Definition: Cubo.cpp:849
void SetCounters()
Set all the counters to zero.
Definition: Cubo.cpp:1088
int SetCoorNumb(double *Pos, int p)
Set the coordination umber to the part p.
Definition: Cubo.cpp:271
void PrintCell(const int c)
Print the content of one cell.
Definition: Cubo.cpp:714
void IncrCurr(int c)
Increase the counters.
Definition: Cubo.cpp:691
void IncrCurrList(const int c)
Increment the pointers.
Definition: Cubo.cpp:835
int IfItCell(int c)
End the iteration in the cell.
Definition: Cubo.cpp:687
The previuos and consecutive particle in the cell linked list.
Definition: Cubo.h:126
int IfItCell(const int c)
Stop the loop and set the counter to zero.
Definition: Cubo.cpp:427
void AddPart(int p, double *Pos)
Add a part to the cell.
Definition: Cubo.cpp:755
void PrintCell(const int c)
Print the content of a cell.
double CutOff
Cell CutOff.
Definition: Cubo.h:65
void IncrCurr(int c)
Increase the counters.
Definition: Cubo.cpp:819
void Add(int v, int t, double *Pos)
Add the triangle t at the vertex v.
Definition: Cubo.cpp:1036
void Couple(const int c, int *p1, int *p2)
Returns the pointed couple of particle.
Definition: Cubo.cpp:697
DomCell & operator++()
Increase the iterator, deactivated.
Definition: Cubo.cpp:611
int pNPart()
part in the class
Definition: Cubo.cpp:51
void IncrCurrList(const int c)
Increment the current iterators in the cell.
Definition: Cubo.cpp:582
void IncrCurr(int v)
Increment the counter.
Definition: Cubo.cpp:1106
void IncrCurr(const int c)
Increment the current part in the cell.
Definition: Cubo.cpp:435
int pCella(const double Pos[3])
(x,y,z)->n
Definition: Cubo.cpp:192
int Mod10[3]
Module 10 (nx,ny,nz) -> nc.
Definition: Cubo.h:55
double Edge[3]
Box sizes.
Definition: Cubo.h:61
int Next
Next particle in the list.
Definition: Cubo.h:130
void PosVertex(int v, double *Pos)
Position of the vertex v.
int Cell
Which cell it&#39;s belonging.
Definition: Cubo.h:136
void RemPart(int p, double *Pos)
Remove a part from the cell.
Definition: Cubo.cpp:659
int operator[](int col)
Returns a entry.
Definition: Cubo.cpp:627
void Rem(int v, int t)
Remove the triangle t from the vertex v.
DomDecBasics()
Allocate.
Definition: Cubo.cpp:17
int ItCell(const int c)
Iterate in the cell.
Definition: Cubo.cpp:421
Definition: Cubo.h:469
int SwapPart(int p1, double *Pos1, int p2, double *Pos2)
Swap two particles.
Definition: Cubo.cpp:342
void SetCounters(int c)
Set the counters to the initial pointer.
Definition: Cubo.cpp:680
int GetCoorNumb(double *Pos)
Get the coordination number.
Definition: Cubo.cpp:212
int NAllocP
Number of allocated part.
Definition: Cubo.h:69
void SwapPart(int p1, double *Pos1, int p2, double *Pos2)
Swap two parts.
Definition: Cubo.cpp:671
void MovePart(int p, double *OldPos, double *NewPos)
Move a part.
Definition: Cubo.cpp:664
For every cube how many particles and which are the first and the last.
Definition: Cubo.h:142
double InvEdge[3]
Inverse box size.
Definition: Cubo.h:63
void Erase()
Erase the content of the cells.
Definition: Cubo.cpp:649
int pNCell()
Print the number of cells.
Definition: Cubo.cpp:55
void MovePart(const int p, double *OldPos, double *NewPos)
Move a particle form the cell c1 to the cell c2.
Definition: Cubo.cpp:403
int IfItCouple(const int c)
Tell when the loop is over.
Definition: Cubo.cpp:830
int Prev
Previous particle in the list.
Definition: Cubo.h:132
The previuos and consecutive particle in the cell linked list.
Definition: Cubo.h:114
int GetVertex(double *Pos)
Correspondent vertes for Pos.
Definition: Cubo.cpp:1025
int ItCell(int c)
Value of the current particle in the cell c.
Definition: Cubo.cpp:822
void PrintCells()
Print the content of every cell.
Definition: Cubo.cpp:720
void Couple(const int c, int *p1, int *p2)
Iterates along all couples.
Definition: Cubo.cpp:565
int GetCell(double *Pos, int *NeiList)
Get the number of the neighbouring cells for a ghost particle.
Definition: Cubo.cpp:59