41 {
42 int mirrorDirection0;
43 int mirrorDirection1;
44 int mirrorDirection2;
45 int mult = 2 / (discreteNormalX*discreteNormalX + discreteNormalY*discreteNormalY + discreteNormalZ*discreteNormalZ);
46 int reflectionPop[DESCRIPTOR::q] = {0};
47 for (int iPop = 1; iPop < DESCRIPTOR::q; iPop++) {
48 reflectionPop[iPop] = 0;
49
50 int scalarProduct = descriptors::c<DESCRIPTOR>(iPop,0)*discreteNormalX + descriptors::c<DESCRIPTOR>(iPop,1)*discreteNormalY + descriptors::c<DESCRIPTOR>(iPop,2)*discreteNormalZ;
51 if (scalarProduct < 0) {
52
53 if (mult == 0) {
54 mirrorDirection0 = -descriptors::c<DESCRIPTOR>(iPop,0);
55 mirrorDirection1 = -descriptors::c<DESCRIPTOR>(iPop,1);
56 mirrorDirection2 = -descriptors::c<DESCRIPTOR>(iPop,2);
57 }
58 else {
59 mirrorDirection0 = descriptors::c<DESCRIPTOR>(iPop,0) - mult*
scalarProduct*discreteNormalX;
60 mirrorDirection1 = descriptors::c<DESCRIPTOR>(iPop,1) - mult*
scalarProduct*discreteNormalY;
61 mirrorDirection2 = descriptors::c<DESCRIPTOR>(iPop,2) - mult*
scalarProduct*discreteNormalZ;
62 }
63
64 for (int i = 1; i < DESCRIPTOR::q; i++) {
65 if (descriptors::c<DESCRIPTOR>(i,0)==mirrorDirection0
66 && descriptors::c<DESCRIPTOR>(i,1)==mirrorDirection1
67 && descriptors::c<DESCRIPTOR>(i,2)==mirrorDirection2) {
68 reflectionPop[iPop] = i;
69 break;
70 }
71 }
72 }
73 }
74 for (int iPop = 1; iPop < DESCRIPTOR::q; iPop++) {
75 if (reflectionPop[iPop]!=0) {
76 cell[iPop] = cell[reflectionPop[iPop]];
77 }
78 }
79 }
T scalarProduct(const T u1[d], const T u2[d])