| 124 | 
  | 
                Vector3d Rij = beads[i].pos - beads[j].pos; | 
| 125 | 
  | 
                double rij = Rij.length(); | 
| 126 | 
  | 
                double rij2 = rij * rij; | 
| 127 | 
< | 
                double sumSigma2OverRij2 = ((beads[i].radius*beads[i].radius) + (beads[i].radius*beads[i].radius)) / rij2;                 | 
| 127 | 
> | 
                double sumSigma2OverRij2 = ((beads[i].radius*beads[i].radius) + (beads[j].radius*beads[j].radius)) / rij2;                 | 
| 128 | 
  | 
                Mat3x3d tmpMat; | 
| 129 | 
  | 
                tmpMat = outProduct(Rij, Rij) / rij2; | 
| 130 | 
  | 
                double constant = 8.0 * NumericConstant::PI * viscosity * rij; | 
| 141 | 
  | 
 | 
| 142 | 
  | 
    //invert B Matrix | 
| 143 | 
  | 
    invertMatrix(B, C); | 
| 144 | 
< | 
 | 
| 144 | 
> | 
     | 
| 145 | 
  | 
    //prepare U Matrix relative to arbitrary origin O(0.0, 0.0, 0.0) | 
| 146 | 
  | 
    std::vector<Mat3x3d> U; | 
| 147 | 
  | 
    for (int i = 0; i < nbeads; ++i) { | 
| 287 | 
  | 
                Vector3d Rij = beads[i].pos - beads[j].pos; | 
| 288 | 
  | 
                double rij = Rij.length(); | 
| 289 | 
  | 
                double rij2 = rij * rij; | 
| 290 | 
< | 
                double sumSigma2OverRij2 = ((beads[i].radius*beads[i].radius) + (beads[i].radius*beads[i].radius)) / rij2;                 | 
| 290 | 
> | 
                double sumSigma2OverRij2 = ((beads[i].radius*beads[i].radius) + (beads[j].radius*beads[j].radius)) / rij2;                 | 
| 291 | 
  | 
                Mat3x3d tmpMat; | 
| 292 | 
  | 
                tmpMat = outProduct(Rij, Rij) / rij2; | 
| 293 | 
  | 
                double constant = 8.0 * NumericConstant::PI * viscosity * rij; |