Line data    Source code

       1             : #include "Compartment.hpp"
       2             : #include "CollisionDetector.hpp"
       3             : #include "CollisionHandler.hpp"
       4             : #include "Disc.hpp"
       5             : #include "MathUtils.hpp"
       6             : #include "ReactionEngine.hpp"
       7             : 
       8             : namespace cell
       9             : {
      10             : 
      11          21 : Compartment::Compartment(Compartment* parent, Membrane membrane, SimulationContext simulationContext)
      12          21 :     : parent_(parent)
      13          21 :     , membrane_(std::move(membrane))
      14          21 :     , simulationContext_(std::move(simulationContext))
      15          42 :     , collisionDetector_(simulationContext_.discTypeRegistry, simulationContext_.membraneTypeRegistry)
      16             : {
      17          21 :     membrane_.setCompartment(this);
      18          21 :     collisionDetector_.setParams(CollisionDetector::Params{.discs = &discs_,
      19          21 :                                                            .membranes = &membranes_,
      20          21 :                                                            .intrudingDiscs = &intrudingDiscs_,
      21          21 :                                                            .containingMembrane = &membrane_});
      22          21 : }
      23             : 
      24          21 : Compartment::~Compartment() = default;
      25             : 
      26      362942 : const Membrane& Compartment::getMembrane() const
      27             : {
      28      362942 :     return membrane_;
      29             : }
      30             : 
      31           0 : void Compartment::setDiscs(std::vector<Disc>&& discs)
      32             : {
      33           0 :     discs_ = std::move(discs);
      34           0 : }
      35             : 
      36         323 : void Compartment::addDisc(Disc disc)
      37             : {
      38         323 :     discs_.push_back(std::move(disc));
      39         323 : }
      40             : 
      41          38 : const std::vector<Disc>& Compartment::getDiscs() const
      42             : {
      43          38 :     return discs_;
      44             : }
      45             : 
      46           5 : void Compartment::addIntrudingDisc(Disc* disc, const Compartment* source, bool shouldBeCaptured)
      47             : {
      48           5 :     intrudingDiscs_.push_back(disc);
      49           5 :     intruderCaptureStatus_.push_back(static_cast<char>(shouldBeCaptured));
      50             : 
      51           5 :     if (compartments_.size() <= 1)
      52           5 :         return;
      53             : 
      54           0 :     const auto discRadius = simulationContext_.discTypeRegistry.getByID(disc->getTypeID()).getRadius();
      55           0 :     for (auto& compartment : compartments_)
      56             :     {
      57           0 :         if (compartment.get() == source)
      58           0 :             continue;
      59             : 
      60             :         // TODO In case of many compartments, make search more efficient
      61             :         const auto* membraneType =
      62           0 :             &simulationContext_.membraneTypeRegistry.getByID(compartment->getMembrane().getTypeID());
      63           0 :         if (mathutils::circlesOverlap(disc->getPosition(), discRadius, compartment->getMembrane().getPosition(),
      64             :                                       membraneType->getRadius()))
      65           0 :             compartment->addIntrudingDisc(disc, source, shouldBeCaptured);
      66             :     }
      67             : }
      68             : 
      69      462366 : std::vector<std::unique_ptr<Compartment>>& Compartment::getCompartments()
      70             : {
      71      462366 :     return compartments_;
      72             : }
      73             : 
      74          14 : const std::vector<std::unique_ptr<Compartment>>& Compartment::getCompartments() const
      75             : {
      76          14 :     return compartments_;
      77             : }
      78             : 
      79           3 : const Compartment* Compartment::getParent() const
      80             : {
      81           3 :     return parent_;
      82             : }
      83             : 
      84           9 : void Compartment::update(double dt)
      85             : {
      86             :     // TODO Remove recursing twice and just accept destroyed discs at the end of update?
      87           9 :     bimolecularUpdate();
      88           9 :     unimolecularUpdate(dt);
      89           9 : }
      90             : 
      91          13 : void Compartment::bimolecularUpdate()
      92             : {
      93          13 :     allocateMemoryForIntruders();
      94          13 :     auto discMembraneCollisions = detectDiscMembraneCollisions();
      95          13 :     registerIntruders(discMembraneCollisions);
      96             : 
      97          17 :     for (auto& compartment : compartments_)
      98           4 :         compartment->bimolecularUpdate();
      99             : 
     100          13 :     auto discDiscCollisions = detectDiscDiscCollisions();
     101          13 :     simulationContext_.collisionHandler.resolveCollisions(discMembraneCollisions);
     102          13 :     simulationContext_.collisionHandler.resolveCollisions(discDiscCollisions);
     103          13 :     simulationContext_.reactionEngine.applyBimolecularReactions(discDiscCollisions, newDiscs_);
     104             : 
     105          13 :     captureIntruders();
     106          13 : }
     107             : 
     108          13 : void Compartment::unimolecularUpdate(double dt)
     109             : {
     110          17 :     for (auto& compartment : compartments_)
     111           4 :         compartment->unimolecularUpdate(dt);
     112             : 
     113          13 :     moveDiscsAndCleanUp(dt);
     114          13 : }
     115             : 
     116          13 : void Compartment::allocateMemoryForIntruders()
     117             : {
     118          13 :     if (intruderAllocationCount_ == 0)
     119          13 :         return;
     120             : 
     121           0 :     discs_.reserve(
     122           0 :         std::max(2 * discs_.capacity(), discs_.capacity() + intruderAllocationCount_ * 2)); // Let's be generous here
     123           0 :     intruderAllocationCount_ = 0;
     124             : }
     125             : 
     126          11 : Compartment* Compartment::createSubCompartment(Membrane membrane)
     127             : {
     128          11 :     auto compartment = std::make_unique<Compartment>(this, std::move(membrane), simulationContext_);
     129          11 :     membranes_.push_back(compartment->getMembrane());
     130          11 :     compartments_.push_back(std::move(compartment));
     131          11 :     collisionDetector_.buildMembraneIndex();
     132             : 
     133          22 :     return compartments_.back().get();
     134          11 : }
     135             : 
     136          13 : std::vector<cell::CollisionDetector::Collision> Compartment::detectDiscMembraneCollisions()
     137             : {
     138          13 :     collisionDetector_.buildDiscIndex();
     139          13 :     return collisionDetector_.detectDiscMembraneCollisions();
     140             : }
     141             : 
     142          13 : std::vector<cell::CollisionDetector::Collision> Compartment::detectDiscDiscCollisions()
     143             : {
     144          13 :     collisionDetector_.addIntrudingDiscsToIndex();
     145          13 :     auto collisions = collisionDetector_.detectDiscDiscCollisions();
     146             : 
     147          13 :     return collisions;
     148             : }
     149             : 
     150          13 : void Compartment::registerIntruders(const std::vector<CollisionDetector::Collision>& discMembraneCollisions)
     151             : {
     152             :     // Assumptions:
     153             :     // - Called directly after collision detection, so no destroyed discs yet
     154             :     // - parent_ can't be nullptr because the outermost membrane shouldn't be permeable for anything so
     155             :     // collision.allowedToPass would always be false
     156             :     // - Also this function only expects DiscContainingMembrane and DiscChildMembrane collisions
     157             : 
     158          21 :     for (const auto& collision : discMembraneCollisions)
     159             :     {
     160           8 :         if (!collision.allowedToPass)
     161           3 :             continue;
     162             : 
     163           5 :         const auto discRadius = simulationContext_.discTypeRegistry.getByID(collision.disc->getTypeID()).getRadius();
     164             : 
     165           5 :         const auto& membranePosition = collision.membrane->getPosition();
     166             :         const auto membraneRadius =
     167           5 :             simulationContext_.membraneTypeRegistry.getByID(collision.membrane->getTypeID()).getRadius();
     168             : 
     169           5 :         if (collision.type == CollisionDetector::CollisionType::DiscContainingMembrane)
     170             :         {
     171             :             // collision.membrane is the membrane of this compartment
     172             :             const bool shouldBeCaptured =
     173           1 :                 !mathutils::circlesOverlap(collision.disc->getPosition(), discRadius, membranePosition, membraneRadius);
     174             : 
     175           1 :             parent_->addIntrudingDisc(collision.disc, this, shouldBeCaptured);
     176             :         }
     177             :         else
     178             :         {
     179             :             // collision.membrane is the child membrane the disc collided with
     180           4 :             const bool shouldBeCaptured = mathutils::circleIsFullyContainedByCircle(
     181           4 :                 collision.disc->getPosition(), discRadius, membranePosition, membraneRadius);
     182             : 
     183           4 :             collision.membrane->getCompartment()->addIntrudingDisc(collision.disc, nullptr, shouldBeCaptured);
     184             :         }
     185             :     }
     186          13 : }
     187             : 
     188          13 : void Compartment::captureIntruders()
     189             : {
     190             :     // Intruders can safely be added without reallocation (which would invalidate any pointers/intruders in other
     191             :     // compartments)
     192          13 :     if (discs_.capacity() >= discs_.size() + intrudingDiscs_.size())
     193             :     {
     194           9 :         for (std::size_t i = 0; i < intrudingDiscs_.size(); ++i)
     195             :         {
     196           0 :             auto& intruder = intrudingDiscs_[i];
     197           0 :             if (intruder->isMarkedDestroyed())
     198           0 :                 continue;
     199             : 
     200           0 :             if (intruderCaptureStatus_[i])
     201             :             {
     202           0 :                 discs_.push_back(*intruder);
     203           0 :                 intruder->markDestroyed();
     204             :             }
     205             :         }
     206             :     }
     207             :     else
     208           4 :         intruderAllocationCount_ = intrudingDiscs_.size();
     209             : 
     210          13 :     intrudingDiscs_.clear();
     211          13 :     intruderCaptureStatus_.clear();
     212          13 : }
     213             : 
     214          13 : void Compartment::moveDiscsAndCleanUp(double dt)
     215             : {
     216          39 :     for (std::size_t i = 0; i < discs_.size(); ++i)
     217             :     {
     218          26 :         auto& disc = discs_[i];
     219          26 :         simulationContext_.reactionEngine.applyUnimolecularReactions(disc, dt, newDiscs_);
     220             : 
     221          26 :         if (disc.isMarkedDestroyed())
     222             :         {
     223          14 :             discs_[i] = std::move(discs_.back());
     224          14 :             discs_.pop_back();
     225          14 :             --i;
     226          14 :             continue;
     227             :         }
     228             : 
     229          12 :         disc.move(disc.getVelocity() * dt);
     230             :     }
     231             : 
     232          23 :     for (auto& disc : newDiscs_)
     233          10 :         disc.move(disc.getVelocity() * dt);
     234             : 
     235          13 :     if (!newDiscs_.empty())
     236             :     {
     237           6 :         discs_.insert(discs_.end(), newDiscs_.begin(), newDiscs_.end());
     238           6 :         newDiscs_.clear();
     239             :     }
     240          13 : }
     241             : 
     242             : } // namespace cell