Line data    Source code

       1             : #include "ReactionEngine.hpp"
       2             : #include "Disc.hpp"
       3             : #include "MathUtils.hpp"
       4             : #include "Reaction.hpp"
       5             : #include "ReactionTable.hpp"
       6             : 
       7             : namespace cell
       8             : {
       9             : 
      10          11 : ReactionEngine::ReactionEngine(const DiscTypeRegistry& discTypeRegistry, const ReactionTable& reactionTable)
      11          11 :     : discTypeRegistry_(discTypeRegistry)
      12             : {
      13          11 :     combineReactionsIntoSingleMaps(reactionTable);
      14          11 : }
      15             : 
      16           1 : Disc ReactionEngine::transformationReaction(Disc* educt, DiscTypeID productID) const
      17             : {
      18           1 :     Disc product(*educt);
      19           1 :     product.setType(productID);
      20           1 :     educt->markDestroyed();
      21             : 
      22           1 :     return product;
      23             : }
      24             : 
      25           1 : std::pair<Disc, Disc> ReactionEngine::decompositionReaction(Disc* educt, DiscTypeID product1ID,
      26             :                                                             DiscTypeID product2ID) const
      27             : {
      28           1 :     double v = mathutils::abs(educt->getVelocity());
      29           1 :     if (v == 0)
      30             :     {
      31             :         // If the disc is stationary and wants to split apart, we'll give it a random velocity to do so
      32           0 :         educt->setVelocity(
      33           0 :             Vector2d{mathutils::getRandomNumber<double>(-10, 10), mathutils::getRandomNumber<double>(-10, 10)});
      34           0 :         v = mathutils::abs(educt->getVelocity());
      35             :     }
      36             : 
      37           1 :     const Vector2d eductNormalizedVelocity = educt->getVelocity() / v;
      38           1 :     const Vector2d n{-eductNormalizedVelocity.y, eductNormalizedVelocity.x};
      39             : 
      40           1 :     Disc product1(*educt);
      41           1 :     Disc product2(product2ID);
      42             : 
      43           1 :     product1.setType(product1ID);
      44           1 :     product1.setVelocity(v * n);
      45             : 
      46           1 :     product2.setPosition(educt->getPosition());
      47           1 :     product2.setVelocity(-v * n);
      48             : 
      49           1 :     const auto R1 = discTypeRegistry_.getByID(product1.getTypeID()).getRadius();
      50           1 :     const auto R2 = discTypeRegistry_.getByID(product2.getTypeID()).getRadius();
      51           1 :     const auto overlap = R1 + R2 + 1e-6; // Discs at same position always have maximum overlap R1 + R2
      52             : 
      53           1 :     product1.move(0.5 * overlap * n);
      54           1 :     product2.move(-0.5 * overlap * n);
      55             : 
      56           1 :     educt->markDestroyed();
      57             : 
      58           2 :     return std::make_pair(std::move(product1), std::move(product2));
      59             : }
      60             : 
      61           5 : Disc ReactionEngine::combinationReaction(Disc* educt1, Disc* educt2, DiscTypeID productID) const
      62             : {
      63           5 :     const auto& resultType = discTypeRegistry_.getByID(productID);
      64           5 :     const auto& d1Type = discTypeRegistry_.getByID(educt1->getTypeID());
      65           5 :     const auto& d2Type = discTypeRegistry_.getByID(educt2->getTypeID());
      66             : 
      67           5 :     Disc newDisc(productID);
      68          10 :     newDisc.setVelocity((d1Type.getMass() * educt1->getVelocity() + d2Type.getMass() * educt2->getVelocity()) /
      69          10 :                         resultType.getMass());
      70           5 :     newDisc.setPosition((educt1->getPosition() + educt2->getPosition()) / 2.0);
      71             : 
      72           5 :     educt1->markDestroyed();
      73           5 :     educt2->markDestroyed();
      74             : 
      75           5 :     return newDisc;
      76             : }
      77             : 
      78           1 : std::pair<Disc, Disc> ReactionEngine::exchangeReaction(Disc* educt1, Disc* educt2, DiscTypeID product1ID,
      79             :                                                        DiscTypeID product2ID) const
      80             : {
      81           1 :     const auto* d1Type = &discTypeRegistry_.getByID(educt1->getTypeID());
      82           1 :     const auto* d2Type = &discTypeRegistry_.getByID(educt2->getTypeID());
      83           1 :     const auto* product1Type = &discTypeRegistry_.getByID(product1ID);
      84           1 :     const auto* product2Type = &discTypeRegistry_.getByID(product2ID);
      85             : 
      86             :     // Sort both product types and educt discs by radius
      87             :     // Now the smallest/largest disc gets the smallest/largest product type
      88             : 
      89           1 :     if (product1Type->getRadius() > product2Type->getRadius())
      90             :     {
      91           0 :         std::swap(product1Type, product2Type);
      92           0 :         std::swap(product1ID, product2ID);
      93             :     }
      94           1 :     if (d1Type->getRadius() > d2Type->getRadius())
      95             :     {
      96           0 :         std::swap(educt1, educt2);
      97           0 :         std::swap(d1Type, d2Type);
      98             :     }
      99             : 
     100             :     // Prefer the assignment that keeps
     101             :     // as many discs as possible with their original type.
     102             : 
     103           1 :     int leaveAsIs = (educt1->getTypeID() == product1ID) + (educt2->getTypeID() == product2ID);
     104           1 :     int swapAgain = (educt1->getTypeID() == product2ID) + (educt2->getTypeID() == product1ID);
     105             : 
     106           1 :     if (swapAgain > leaveAsIs)
     107             :     {
     108           0 :         std::swap(product1Type, product2Type);
     109           0 :         std::swap(product1ID, product2ID);
     110             :     }
     111             : 
     112           1 :     Disc product1(*educt1);
     113           1 :     Disc product2(*educt2);
     114             : 
     115           1 :     product1.scaleVelocity(std::sqrt(d1Type->getMass() / product1Type->getMass()));
     116           1 :     product1.setType(product1ID);
     117             : 
     118           1 :     product2.scaleVelocity(std::sqrt(d2Type->getMass() / product2Type->getMass()));
     119           1 :     product2.setType(product2ID);
     120             : 
     121           1 :     educt1->markDestroyed();
     122           1 :     educt2->markDestroyed();
     123             : 
     124           2 :     return std::make_pair(std::move(product1), std::move(product2));
     125             : }
     126             : 
     127          26 : void ReactionEngine::applyUnimolecularReactions(Disc& disc, double dt, std::vector<Disc>& newDiscs) const
     128             : {
     129          26 :     const Reaction* reaction = selectUnimolecularReaction(disc.getTypeID(), dt);
     130          26 :     if (!reaction)
     131          24 :         return;
     132             : 
     133           2 :     if (reaction->getType() == Reaction::Type::Transformation)
     134             :     {
     135           1 :         auto product = transformationReaction(&disc, reaction->getProduct1());
     136           1 :         newDiscs.push_back(std::move(product));
     137             :     }
     138             :     else
     139             :     {
     140           1 :         auto products = decompositionReaction(&disc, reaction->getProduct1(), reaction->getProduct2());
     141           1 :         newDiscs.push_back(std::move(products.first));
     142           1 :         newDiscs.push_back(std::move(products.second));
     143             :     }
     144             : }
     145             : 
     146          13 : void ReactionEngine::applyBimolecularReactions(const std::vector<CollisionDetector::Collision>& collisions,
     147             :                                                std::vector<Disc>& newDiscs) const
     148             : {
     149          20 :     for (const auto& collision : collisions)
     150             :     {
     151           7 :         if (collision.invalidatedByDestroyedDiscs())
     152           0 :             continue;
     153             : 
     154             :         const Reaction* reaction =
     155           7 :             selectBimolecularReaction(std::minmax(collision.disc->getTypeID(), collision.otherDisc->getTypeID()));
     156           7 :         if (!reaction)
     157           1 :             continue;
     158             : 
     159           6 :         if (reaction->getType() == Reaction::Type::Combination)
     160             :         {
     161           5 :             auto product = combinationReaction(collision.disc, collision.otherDisc, reaction->getProduct1());
     162           5 :             newDiscs.push_back(std::move(product));
     163             :         }
     164             :         else
     165             :         {
     166             :             auto products =
     167           1 :                 exchangeReaction(collision.disc, collision.otherDisc, reaction->getProduct1(), reaction->getProduct2());
     168           1 :             newDiscs.push_back(std::move(products.first));
     169           1 :             newDiscs.push_back(std::move(products.second));
     170             :         }
     171             :     }
     172          13 : }
     173             : 
     174          26 : const Reaction* ReactionEngine::selectUnimolecularReaction(const DiscTypeID& key, double dt) const
     175             : {
     176          26 :     return selectReaction(
     177          26 :         unimolecularReactions_, key, [&](const Reaction& reaction)
     178          54 :         { return mathutils::getRandomNumber<double>(0, 1) <= 1 - std::pow(1 - reaction.getProbability(), dt); });
     179             : }
     180             : 
     181           7 : const Reaction* ReactionEngine::selectBimolecularReaction(const std::pair<DiscTypeID, DiscTypeID>& key) const
     182             : {
     183          14 :     return selectReaction(bimolecularReactions_, key, [](const Reaction& reaction)
     184          20 :                           { return mathutils::getRandomNumber<double>(0, 1) <= reaction.getProbability(); });
     185             : }
     186             : 
     187          11 : void ReactionEngine::combineReactionsIntoSingleMaps(const ReactionTable& reactionTable)
     188             : {
     189          55 :     for (const auto& table : {reactionTable.getTransformations(), reactionTable.getDecompositions()})
     190             :     {
     191          24 :         for (const auto& [educt, reactions] : table)
     192           2 :             unimolecularReactions_[educt].insert(unimolecularReactions_[educt].end(), reactions.begin(),
     193             :                                                  reactions.end());
     194          33 :     }
     195             : 
     196          55 :     for (const auto& table : {reactionTable.getCombinations(), reactionTable.getExchanges()})
     197             :     {
     198          27 :         for (const auto& [educts, reactions] : table)
     199           5 :             bimolecularReactions_[educts].insert(bimolecularReactions_[educts].end(), reactions.begin(),
     200             :                                                  reactions.end());
     201          33 :     }
     202             : 
     203             :     // Random shuffle all reactions to avoid a bias
     204          11 :     std::minstd_rand rng{std::random_device{}()};
     205             : 
     206          13 :     for (auto& [educt, reactions] : unimolecularReactions_)
     207           2 :         std::shuffle(reactions.begin(), reactions.end(), rng);
     208             : 
     209          16 :     for (auto& [educts, reactions] : bimolecularReactions_)
     210           5 :         std::shuffle(reactions.begin(), reactions.end(), rng);
     211          11 : }
     212             : 
     213             : } // namespace cell