(1) // Scalar interaction |
(2) const RealType dx = positionsX[idxTarget] − positionsX[idxSource]; |
(3) const RealType dy = positionsY[idxTarget] − positionsY[idxSource]; |
(4) const RealType dz = positionsZ[idxTarget] − positionsZ[idxSource]; |
(5) |
(6) const RealType distance = std::sqrt(dxdx + dydy + dzdz); |
(7) const RealType inv_distance = 1/distance; |
(8) |
(9) if(distance < cutDistance) |
(10) potentials[idxTarget] += ( inv_distance physicalValues[idxSource] ); |
(11) potentials[idxSource] += ( inv_distance physicalValues[idxTarget] ); |
(12) |
(13) else |
(14) potentials[idxTarget] += ( inv_distance (physicalValues[idxSource]−constantIfCut) ); |
(15) potentials[idxSource] += ( inv_distance (physicalValues[idxTarget]−constantIfCut) ); |
(16) |
(17) |
(18) // Vectorized using advanced branch manager |
(19) const VecType dx = targetX − VecType(&positionsX[idxSource]); |
(20) const VecType dy = targetY − VecType(&positionsY[idxSource]); |
(21) const VecType dz = targetZ − VecType(&positionsZ[idxSource]); |
(22) |
(23) const VecType distance = VecType(dxdx + dydy + dzdz).sqrt(); |
(24) const VecType inv_distance = VecOne/distance; |
(25) |
(26) const typename VecType::MaskType testRes = (distance < VecCutDistance); |
(27) |
(28) const VecType sourcesPhysicalValue = VecType(&physicalValues[idxSource]); |
(29) |
(30) targetPotential += inv_distance VecType::If(testRes).Then([&]() |
(31) return sourcesPhysicalValue; |
(32) ).Else([&]() |
(33) return sourcesPhysicalValue−VecConstantIfCut; |
(34) ); |
(35) const VecType resSource = inv_distance VecType::If(testRes).Then([&]() |
(36) return targetPhysicalValue; |
(37) ).Else([&]() |
(38) return targetPhysicalValue−VecConstantIfCut; |
(39) ); |
(40) |
(41) const VecType currentSource = VecType(&potentials[idxSource]); |
(42) (resSource+currentSource).storeInArray(&potentials[idxSource]); |