45 {
46 using DESCRIPTOR = typename CELL::descriptor_t;
47 int reflectionPop[DESCRIPTOR::q];
48 int mirrorDirection0;
49 int mirrorDirection1;
50 int mirrorDirection2;
51 int mult = 2 / (NX*NX + NY*NY + NZ*NZ);
52 reflectionPop[0] =0;
53 for (int iPop = 1; iPop < DESCRIPTOR::q; iPop++) {
54 reflectionPop[iPop] = 0;
55
56 int scalarProduct = descriptors::c<DESCRIPTOR>(iPop,0)*NX + descriptors::c<DESCRIPTOR>(iPop,1)*NY + descriptors::c<DESCRIPTOR>(iPop,2)*NZ;
57 if ( scalarProduct < 0) {
58
59 if (mult == 0) {
60 mirrorDirection0 = -descriptors::c<DESCRIPTOR>(iPop,0);
61 mirrorDirection1 = -descriptors::c<DESCRIPTOR>(iPop,1);
62 mirrorDirection2 = -descriptors::c<DESCRIPTOR>(iPop,2);
63 }
64 else {
65 mirrorDirection0 = descriptors::c<DESCRIPTOR>(iPop,0) - mult*
scalarProduct*NX;
66 mirrorDirection1 = descriptors::c<DESCRIPTOR>(iPop,1) - mult*
scalarProduct*NY;
67 mirrorDirection2 = descriptors::c<DESCRIPTOR>(iPop,2) - mult*
scalarProduct*NZ;
68 }
69
70
71 for (int i = 1; i < DESCRIPTOR::q; i++) {
72 if (descriptors::c<DESCRIPTOR>(i,0)==mirrorDirection0
73 && descriptors::c<DESCRIPTOR>(i,1)==mirrorDirection1
74 && descriptors::c<DESCRIPTOR>(i,2)==mirrorDirection2) {
75 reflectionPop[iPop] = i;
76 break;
77 }
78 }
79 }
80 }
81 for (int iPop = 1; iPop < DESCRIPTOR::q ; ++iPop) {
82 if (reflectionPop[iPop]!=0) {
83
84 x_b[iPop] = x_b[reflectionPop[iPop]];
85 }
86 }
87 }
T scalarProduct(const T u1[d], const T u2[d])