From c19e0ad7bf7be96aa7ce23d0a6da7e703cd2f4fb Mon Sep 17 00:00:00 2001
From: Sebastian Eibl <sebastian.eibl@fau.de>
Date: Tue, 26 May 2020 11:43:58 +0200
Subject: [PATCH] [UPDATE] thermal extension fully supports linear expansion

---
 .../HeatConduction/ContactDetection.h         | 78 ++-----------------
 .../HeatConduction/HeatConduction.cpp         | 15 +++-
 .../HeatConduction/ThermalExpansion.h         | 43 ++++++++--
 python/mesa_pd.py                             |  2 +-
 src/mesa_pd/data/ParticleAccessor.h           | 14 ++--
 src/mesa_pd/data/ParticleStorage.h            | 44 +++++------
 src/mesa_pd/mpi/notifications/ParseMessage.h  |  2 +-
 .../notifications/ParticleCopyNotification.h  |  8 +-
 .../ParticleGhostCopyNotification.h           |  8 +-
 .../ParticleUpdateNotification.h              |  6 +-
 10 files changed, 100 insertions(+), 120 deletions(-)

diff --git a/apps/showcases/HeatConduction/ContactDetection.h b/apps/showcases/HeatConduction/ContactDetection.h
index f1c0fd5d1..c5ed90565 100644
--- a/apps/showcases/HeatConduction/ContactDetection.h
+++ b/apps/showcases/HeatConduction/ContactDetection.h
@@ -32,30 +32,13 @@
 namespace walberla {
 namespace mesa_pd {
 
-class ContactDetection
+class ContactDetection : public collision_detection::AnalyticContactDetection
 {
 public:
-   size_t& getIdx1() {return idx1_;}
-   size_t& getIdx2() {return idx2_;}
-   Vec3&   getContactPoint() {return contactPoint_;}
-   Vec3&   getContactNormal() {return contactNormal_;}
-   real_t& getPenetrationDepth() {return penetrationDepth_;}
-
-   const size_t& getIdx1() const {return idx1_;}
-   const size_t& getIdx2() const {return idx2_;}
-   const Vec3&   getContactPoint() const {return contactPoint_;}
-   const Vec3&   getContactNormal() const {return contactNormal_;}
-   const real_t& getPenetrationDepth() const {return penetrationDepth_;}
-
-   real_t& getContactThreshold() {return contactThreshold_;}
-
-   template <typename GEO1_T, typename GEO2_T, typename Accessor>
-   bool operator()( const size_t idx1,
-                    const size_t idx2,
-                    const GEO1_T& geo1,
-                    const GEO2_T& geo2,
-                    Accessor& ac);
+   ///import all collision detection functions from AnalyticContactDetection
+   using AnalyticContactDetection::operator();
 
+   ///overwrite functions that should be different
    template <typename Accessor>
    bool operator()( const size_t idx1,
                     const size_t idx2,
@@ -69,34 +52,8 @@ public:
                     const data::Sphere& s,
                     const data::HalfSpace& p,
                     Accessor& ac );
-
-   template <typename Accessor>
-   bool operator()( const size_t idx1,
-                    const size_t idx2,
-                    const data::HalfSpace& p,
-                    const data::Sphere& s,
-                    Accessor& ac);
-
-private:
-   size_t idx1_;
-   size_t idx2_;
-   Vec3   contactPoint_;
-   Vec3   contactNormal_;
-   real_t penetrationDepth_;
-
-   real_t contactThreshold_ = real_t(0.0);
 };
 
-template <typename GEO1_T, typename GEO2_T, typename Accessor>
-inline bool ContactDetection::operator()( const size_t /*idx1*/,
-                                          const size_t /*idx2*/,
-                                          const GEO1_T& /*geo1*/,
-                                          const GEO2_T& /*geo2*/,
-                                          Accessor& /*ac*/)
-{
-   WALBERLA_ABORT("Collision not implemented!")
-}
-
 template <typename Accessor>
 inline bool ContactDetection::operator()( const size_t idx1,
                                           const size_t idx2,
@@ -114,9 +71,9 @@ inline bool ContactDetection::operator()( const size_t idx1,
    getIdx1() = idx1;
    getIdx2() = idx2;
    return detectSphereSphereCollision(ac.getPosition(getIdx1()),
-                                      ac.getRadius(getIdx1()),
+                                      ac.getRadiusAtTemperature(getIdx1()),
                                       ac.getPosition(getIdx2()),
-                                      ac.getRadius(getIdx2()),
+                                      ac.getRadiusAtTemperature(getIdx2()),
                                       getContactPoint(),
                                       getContactNormal(),
                                       getPenetrationDepth(),
@@ -136,7 +93,7 @@ inline bool ContactDetection::operator()( const size_t idx1,
    getIdx1() = idx1;
    getIdx2() = idx2;
    return detectSphereHalfSpaceCollision(ac.getPosition(getIdx1()),
-                                         ac.getRadius(getIdx1()),
+                                         ac.getRadiusAtTemperature(getIdx1()),
                                          ac.getPosition(getIdx2()),
                                          p.getNormal(),
                                          getContactPoint(),
@@ -145,26 +102,5 @@ inline bool ContactDetection::operator()( const size_t idx1,
                                          getContactThreshold());
 }
 
-template <typename Accessor>
-inline bool ContactDetection::operator()( const size_t idx1,
-                                          const size_t idx2,
-                                          const data::HalfSpace& p,
-                                          const data::Sphere& s,
-                                          Accessor& ac)
-{
-   return operator()(idx2, idx1, s, p, ac);
-}
-
-inline
-std::ostream& operator<<( std::ostream& os, const ContactDetection& ac )
-{
-   os << "idx1:               " << ac.getIdx1() << "\n" <<
-         "idx2:               " << ac.getIdx2() << "\n" <<
-         "contact point:      " << ac.getContactPoint() << "\n" <<
-         "contact normal:     " << ac.getContactNormal() << "\n" <<
-         "penetration depth:  " << ac.getPenetrationDepth() << std::endl;
-   return os;
-}
-
 } //namespace mesa_pd
 } //namespace walberla
diff --git a/apps/showcases/HeatConduction/HeatConduction.cpp b/apps/showcases/HeatConduction/HeatConduction.cpp
index 72bf5755d..eea1b2e72 100644
--- a/apps/showcases/HeatConduction/HeatConduction.cpp
+++ b/apps/showcases/HeatConduction/HeatConduction.cpp
@@ -18,6 +18,17 @@
 //
 //======================================================================================================================
 
+//======================================================================================================================
+// This showcase combines a classic spring-dashpot interaction model with heat conduction and linear thermal expansion.
+//
+// First, a set of spheres is dropped onto a plate. The setup is periodic in horizontal directions. The plate as well
+// as the spheres shares the same initial temperature. After the spheres have settled, the plate is heated and cooled
+// down in a sinus wave. Through thermal conduction the heat is transported through the pile of spheres. At the same
+// time a linear thermal expansion model makes the spheres grow and shrink. The parameters are chosen for a visual
+// effect and are far to large for a realistic physical simulation. The evolution of the simulation is written to disk
+// as vtk files which can be visualized with paraview.
+//======================================================================================================================
+
 #include <mesa_pd/vtk/ParticleVtkOutput.h>
 
 #include <mesa_pd/collision_detection/BroadPhase.h>
@@ -179,7 +190,7 @@ int main( int argc, char ** argv )
          p->getOwnerRef()             = walberla::mpi::MPIManager::instance()->rank();
          p->getTypeRef()              = 0;
          p->setTemperature( 273 );
-         p->setRadius(radius);
+         p->setRadiusAtTemperature(radius);
       }
    }
    int64_t numParticles = int64_c(ps->size());
@@ -197,7 +208,7 @@ int main( int argc, char ** argv )
    auto vtkWriter       = walberla::vtk::createVTKOutput_PointData(vtkOutput, "particles", 1, "vtk", "simulation_step", false, false);
    vtkOutput->addOutput<data::SelectParticleOwner>("owner");
    vtkOutput->addOutput<data::SelectParticleTemperature>("temperature");
-   vtkOutput->addOutput<data::SelectParticleRadius>("radius");
+   vtkOutput->addOutput<data::SelectParticleRadiusAtTemperature>("radius");
    vtkOutput->addOutput<data::SelectParticleLinearVelocity>("linVel");
    vtkOutput->setParticleSelector([smallSphere](const data::ParticleStorage::iterator& pIt){ return pIt->getShapeID() == smallSphere;});
    vtkDomainOutput->write();
diff --git a/apps/showcases/HeatConduction/ThermalExpansion.h b/apps/showcases/HeatConduction/ThermalExpansion.h
index e553c2692..c56ee0732 100644
--- a/apps/showcases/HeatConduction/ThermalExpansion.h
+++ b/apps/showcases/HeatConduction/ThermalExpansion.h
@@ -30,20 +30,53 @@ namespace kernel {
 class ThermalExpansion
 {
 public:
-   ThermalExpansion() = default;
+   ThermalExpansion(const uint_t numParticleTypes);
+   ThermalExpansion(const ThermalExpansion& other) = default;
+   ThermalExpansion(ThermalExpansion&& other) = default;
+   ThermalExpansion& operator=(const ThermalExpansion& other) = default;
+   ThermalExpansion& operator=(ThermalExpansion&& other) = default;
+
+   void setLinearExpansionCoefficient(const size_t type, const real_t& val);
+   real_t getLinearExpansionCoefficient(const size_t type) const;
 
    template <typename Accessor>
-   void operator()(const size_t i, Accessor& ac) const;
+   void operator()(const size_t p_idx, Accessor& ac) const;
+private:
+       uint_t numParticleTypes_;
+
+   std::vector<real_t> linearExpansionCoefficient_ {};
 };
 
+ThermalExpansion::ThermalExpansion(const uint_t numParticleTypes)
+{
+   numParticleTypes_ = numParticleTypes;
+
+   linearExpansionCoefficient_.resize(numParticleTypes, real_t(0));
+}
+
+
+inline void ThermalExpansion::setLinearExpansionCoefficient(const size_t type, const real_t& val)
+{
+   WALBERLA_ASSERT_LESS( type, numParticleTypes_ );
+   linearExpansionCoefficient_[type] = val;
+}
+
+
+inline real_t ThermalExpansion::getLinearExpansionCoefficient(const size_t type) const
+{
+   WALBERLA_ASSERT_LESS( type, numParticleTypes_ );
+   return linearExpansionCoefficient_[type];
+}
+
 template <typename Accessor>
-inline void ThermalExpansion::operator()(const size_t idx,
+inline void ThermalExpansion::operator()(const size_t p_idx,
                                          Accessor& ac) const
 {
    static_assert(std::is_base_of<data::IAccessor, Accessor>::value, "please provide a valid accessor");
 
-   ac.setRadius(idx, real_t(0.004)  + (ac.getTemperature(idx)-real_t(273)) * real_t(0.002) * real_t(0.001));
-   ac.setInteractionRadius(idx, ac.getRadius(idx));
+   const auto Tc = ac.getTemperature(p_idx)-real_t(273);
+   ac.setRadiusAtTemperature(p_idx, ac.getRadius273K(p_idx) * (real_t(1) + Tc * getLinearExpansionCoefficient(ac.getType(p_idx))));
+   ac.setInteractionRadius(p_idx, ac.getRadiusAtTemperature(p_idx));
 }
 
 } //namespace kernel
diff --git a/python/mesa_pd.py b/python/mesa_pd.py
index dd53b2f4c..a6d7fa151 100755
--- a/python/mesa_pd.py
+++ b/python/mesa_pd.py
@@ -28,7 +28,7 @@ if __name__ == '__main__':
     ps.add_property("torque", "walberla::mesa_pd::Vec3", defValue="real_t(0)", syncMode="NEVER")
     ps.add_property("oldTorque", "walberla::mesa_pd::Vec3", defValue="real_t(0)", syncMode="ON_OWNERSHIP_CHANGE")
 
-    ps.add_property("radius", "walberla::real_t", defValue="real_t(0)", syncMode="ALWAYS")
+    ps.add_property("radiusAtTemperature", "walberla::real_t", defValue="real_t(0)", syncMode="ALWAYS")
 
     ps.add_include("blockforest/BlockForest.h")
     ps.add_property("currentBlock", "blockforest::BlockID", defValue="", syncMode="NEVER")
diff --git a/src/mesa_pd/data/ParticleAccessor.h b/src/mesa_pd/data/ParticleAccessor.h
index 753546900..a00eb26c3 100644
--- a/src/mesa_pd/data/ParticleAccessor.h
+++ b/src/mesa_pd/data/ParticleAccessor.h
@@ -108,9 +108,9 @@ public:
    walberla::mesa_pd::Vec3& getOldTorqueRef(const size_t p_idx) {return ps_->getOldTorqueRef(p_idx);}
    void setOldTorque(const size_t p_idx, walberla::mesa_pd::Vec3 const & v) { ps_->setOldTorque(p_idx, v);}
    
-   walberla::real_t const & getRadius(const size_t p_idx) const {return ps_->getRadius(p_idx);}
-   walberla::real_t& getRadiusRef(const size_t p_idx) {return ps_->getRadiusRef(p_idx);}
-   void setRadius(const size_t p_idx, walberla::real_t const & v) { ps_->setRadius(p_idx, v);}
+   walberla::real_t const & getRadiusAtTemperature(const size_t p_idx) const {return ps_->getRadiusAtTemperature(p_idx);}
+   walberla::real_t& getRadiusAtTemperatureRef(const size_t p_idx) {return ps_->getRadiusAtTemperatureRef(p_idx);}
+   void setRadiusAtTemperature(const size_t p_idx, walberla::real_t const & v) { ps_->setRadiusAtTemperature(p_idx, v);}
    
    blockforest::BlockID const & getCurrentBlock(const size_t p_idx) const {return ps_->getCurrentBlock(p_idx);}
    blockforest::BlockID& getCurrentBlockRef(const size_t p_idx) {return ps_->getCurrentBlockRef(p_idx);}
@@ -273,9 +273,9 @@ public:
    void setOldTorque(const size_t /*p_idx*/, walberla::mesa_pd::Vec3 const & v) { oldTorque_ = v;}
    walberla::mesa_pd::Vec3& getOldTorqueRef(const size_t /*p_idx*/) {return oldTorque_;}
    
-   walberla::real_t const & getRadius(const size_t /*p_idx*/) const {return radius_;}
-   void setRadius(const size_t /*p_idx*/, walberla::real_t const & v) { radius_ = v;}
-   walberla::real_t& getRadiusRef(const size_t /*p_idx*/) {return radius_;}
+   walberla::real_t const & getRadiusAtTemperature(const size_t /*p_idx*/) const {return radiusAtTemperature_;}
+   void setRadiusAtTemperature(const size_t /*p_idx*/, walberla::real_t const & v) { radiusAtTemperature_ = v;}
+   walberla::real_t& getRadiusAtTemperatureRef(const size_t /*p_idx*/) {return radiusAtTemperature_;}
    
    blockforest::BlockID const & getCurrentBlock(const size_t /*p_idx*/) const {return currentBlock_;}
    void setCurrentBlock(const size_t /*p_idx*/, blockforest::BlockID const & v) { currentBlock_ = v;}
@@ -359,7 +359,7 @@ private:
    walberla::mesa_pd::Vec3 angularVelocity_;
    walberla::mesa_pd::Vec3 torque_;
    walberla::mesa_pd::Vec3 oldTorque_;
-   walberla::real_t radius_;
+   walberla::real_t radiusAtTemperature_;
    blockforest::BlockID currentBlock_;
    uint_t type_;
    int nextParticle_;
diff --git a/src/mesa_pd/data/ParticleStorage.h b/src/mesa_pd/data/ParticleStorage.h
index fd1c33cf8..8e248aca1 100644
--- a/src/mesa_pd/data/ParticleStorage.h
+++ b/src/mesa_pd/data/ParticleStorage.h
@@ -86,7 +86,7 @@ public:
       using angularVelocity_type = walberla::mesa_pd::Vec3;
       using torque_type = walberla::mesa_pd::Vec3;
       using oldTorque_type = walberla::mesa_pd::Vec3;
-      using radius_type = walberla::real_t;
+      using radiusAtTemperature_type = walberla::real_t;
       using currentBlock_type = blockforest::BlockID;
       using type_type = uint_t;
       using nextParticle_type = int;
@@ -163,9 +163,9 @@ public:
       oldTorque_type& getOldTorqueRef() {return storage_.getOldTorqueRef(i_);}
       void setOldTorque(oldTorque_type const & v) { storage_.setOldTorque(i_, v);}
       
-      radius_type const & getRadius() const {return storage_.getRadius(i_);}
-      radius_type& getRadiusRef() {return storage_.getRadiusRef(i_);}
-      void setRadius(radius_type const & v) { storage_.setRadius(i_, v);}
+      radiusAtTemperature_type const & getRadiusAtTemperature() const {return storage_.getRadiusAtTemperature(i_);}
+      radiusAtTemperature_type& getRadiusAtTemperatureRef() {return storage_.getRadiusAtTemperatureRef(i_);}
+      void setRadiusAtTemperature(radiusAtTemperature_type const & v) { storage_.setRadiusAtTemperature(i_, v);}
       
       currentBlock_type const & getCurrentBlock() const {return storage_.getCurrentBlock(i_);}
       currentBlock_type& getCurrentBlockRef() {return storage_.getCurrentBlockRef(i_);}
@@ -298,7 +298,7 @@ public:
    using angularVelocity_type = walberla::mesa_pd::Vec3;
    using torque_type = walberla::mesa_pd::Vec3;
    using oldTorque_type = walberla::mesa_pd::Vec3;
-   using radius_type = walberla::real_t;
+   using radiusAtTemperature_type = walberla::real_t;
    using currentBlock_type = blockforest::BlockID;
    using type_type = uint_t;
    using nextParticle_type = int;
@@ -375,9 +375,9 @@ public:
    oldTorque_type& getOldTorqueRef(const size_t idx) {return oldTorque_[idx];}
    void setOldTorque(const size_t idx, oldTorque_type const & v) { oldTorque_[idx] = v; }
    
-   radius_type const & getRadius(const size_t idx) const {return radius_[idx];}
-   radius_type& getRadiusRef(const size_t idx) {return radius_[idx];}
-   void setRadius(const size_t idx, radius_type const & v) { radius_[idx] = v; }
+   radiusAtTemperature_type const & getRadiusAtTemperature(const size_t idx) const {return radiusAtTemperature_[idx];}
+   radiusAtTemperature_type& getRadiusAtTemperatureRef(const size_t idx) {return radiusAtTemperature_[idx];}
+   void setRadiusAtTemperature(const size_t idx, radiusAtTemperature_type const & v) { radiusAtTemperature_[idx] = v; }
    
    currentBlock_type const & getCurrentBlock(const size_t idx) const {return currentBlock_[idx];}
    currentBlock_type& getCurrentBlockRef(const size_t idx) {return currentBlock_[idx];}
@@ -541,7 +541,7 @@ public:
    std::vector<angularVelocity_type> angularVelocity_ {};
    std::vector<torque_type> torque_ {};
    std::vector<oldTorque_type> oldTorque_ {};
-   std::vector<radius_type> radius_ {};
+   std::vector<radiusAtTemperature_type> radiusAtTemperature_ {};
    std::vector<currentBlock_type> currentBlock_ {};
    std::vector<type_type> type_ {};
    std::vector<nextParticle_type> nextParticle_ {};
@@ -580,7 +580,7 @@ ParticleStorage::Particle& ParticleStorage::Particle::operator=(const ParticleSt
    getAngularVelocityRef() = rhs.getAngularVelocity();
    getTorqueRef() = rhs.getTorque();
    getOldTorqueRef() = rhs.getOldTorque();
-   getRadiusRef() = rhs.getRadius();
+   getRadiusAtTemperatureRef() = rhs.getRadiusAtTemperature();
    getCurrentBlockRef() = rhs.getCurrentBlock();
    getTypeRef() = rhs.getType();
    getNextParticleRef() = rhs.getNextParticle();
@@ -616,7 +616,7 @@ ParticleStorage::Particle& ParticleStorage::Particle::operator=(ParticleStorage:
    getAngularVelocityRef() = std::move(rhs.getAngularVelocityRef());
    getTorqueRef() = std::move(rhs.getTorqueRef());
    getOldTorqueRef() = std::move(rhs.getOldTorqueRef());
-   getRadiusRef() = std::move(rhs.getRadiusRef());
+   getRadiusAtTemperatureRef() = std::move(rhs.getRadiusAtTemperatureRef());
    getCurrentBlockRef() = std::move(rhs.getCurrentBlockRef());
    getTypeRef() = std::move(rhs.getTypeRef());
    getNextParticleRef() = std::move(rhs.getNextParticleRef());
@@ -653,7 +653,7 @@ void swap(ParticleStorage::Particle lhs, ParticleStorage::Particle rhs)
    std::swap(lhs.getAngularVelocityRef(), rhs.getAngularVelocityRef());
    std::swap(lhs.getTorqueRef(), rhs.getTorqueRef());
    std::swap(lhs.getOldTorqueRef(), rhs.getOldTorqueRef());
-   std::swap(lhs.getRadiusRef(), rhs.getRadiusRef());
+   std::swap(lhs.getRadiusAtTemperatureRef(), rhs.getRadiusAtTemperatureRef());
    std::swap(lhs.getCurrentBlockRef(), rhs.getCurrentBlockRef());
    std::swap(lhs.getTypeRef(), rhs.getTypeRef());
    std::swap(lhs.getNextParticleRef(), rhs.getNextParticleRef());
@@ -690,7 +690,7 @@ std::ostream& operator<<( std::ostream& os, const ParticleStorage::Particle& p )
          "angularVelocity     : " << p.getAngularVelocity() << "\n" <<
          "torque              : " << p.getTorque() << "\n" <<
          "oldTorque           : " << p.getOldTorque() << "\n" <<
-         "radius              : " << p.getRadius() << "\n" <<
+         "radiusAtTemperature : " << p.getRadiusAtTemperature() << "\n" <<
          "currentBlock        : " << p.getCurrentBlock() << "\n" <<
          "type                : " << p.getType() << "\n" <<
          "nextParticle        : " << p.getNextParticle() << "\n" <<
@@ -797,7 +797,7 @@ inline ParticleStorage::iterator ParticleStorage::create(const id_t& uid)
    angularVelocity_.emplace_back(real_t(0));
    torque_.emplace_back(real_t(0));
    oldTorque_.emplace_back(real_t(0));
-   radius_.emplace_back(real_t(0));
+   radiusAtTemperature_.emplace_back(real_t(0));
    currentBlock_.emplace_back();
    type_.emplace_back(0);
    nextParticle_.emplace_back(-1);
@@ -859,7 +859,7 @@ inline ParticleStorage::iterator ParticleStorage::erase(iterator& it)
    angularVelocity_.pop_back();
    torque_.pop_back();
    oldTorque_.pop_back();
-   radius_.pop_back();
+   radiusAtTemperature_.pop_back();
    currentBlock_.pop_back();
    type_.pop_back();
    nextParticle_.pop_back();
@@ -908,7 +908,7 @@ inline void ParticleStorage::reserve(const size_t size)
    angularVelocity_.reserve(size);
    torque_.reserve(size);
    oldTorque_.reserve(size);
-   radius_.reserve(size);
+   radiusAtTemperature_.reserve(size);
    currentBlock_.reserve(size);
    type_.reserve(size);
    nextParticle_.reserve(size);
@@ -942,7 +942,7 @@ inline void ParticleStorage::clear()
    angularVelocity_.clear();
    torque_.clear();
    oldTorque_.clear();
-   radius_.clear();
+   radiusAtTemperature_.clear();
    currentBlock_.clear();
    type_.clear();
    nextParticle_.clear();
@@ -977,7 +977,7 @@ inline size_t ParticleStorage::size() const
    //WALBERLA_ASSERT_EQUAL( uid_.size(), angularVelocity.size() );
    //WALBERLA_ASSERT_EQUAL( uid_.size(), torque.size() );
    //WALBERLA_ASSERT_EQUAL( uid_.size(), oldTorque.size() );
-   //WALBERLA_ASSERT_EQUAL( uid_.size(), radius.size() );
+   //WALBERLA_ASSERT_EQUAL( uid_.size(), radiusAtTemperature.size() );
    //WALBERLA_ASSERT_EQUAL( uid_.size(), currentBlock.size() );
    //WALBERLA_ASSERT_EQUAL( uid_.size(), type.size() );
    //WALBERLA_ASSERT_EQUAL( uid_.size(), nextParticle.size() );
@@ -1311,13 +1311,13 @@ public:
    walberla::mesa_pd::Vec3 const & operator()(const data::Particle& p) const {return p.getOldTorque();}
 };
 ///Predicate that selects a certain property from a Particle
-class SelectParticleRadius
+class SelectParticleRadiusAtTemperature
 {
 public:
    using return_type = walberla::real_t;
-   walberla::real_t& operator()(data::Particle& p) const {return p.getRadiusRef();}
-   walberla::real_t& operator()(data::Particle&& p) const {return p.getRadiusRef();}
-   walberla::real_t const & operator()(const data::Particle& p) const {return p.getRadius();}
+   walberla::real_t& operator()(data::Particle& p) const {return p.getRadiusAtTemperatureRef();}
+   walberla::real_t& operator()(data::Particle&& p) const {return p.getRadiusAtTemperatureRef();}
+   walberla::real_t const & operator()(const data::Particle& p) const {return p.getRadiusAtTemperature();}
 };
 ///Predicate that selects a certain property from a Particle
 class SelectParticleCurrentBlock
diff --git a/src/mesa_pd/mpi/notifications/ParseMessage.h b/src/mesa_pd/mpi/notifications/ParseMessage.h
index 02b36290f..caf0e0f89 100644
--- a/src/mesa_pd/mpi/notifications/ParseMessage.h
+++ b/src/mesa_pd/mpi/notifications/ParseMessage.h
@@ -117,7 +117,7 @@ void ParseMessage::operator()(int sender,
       pIt->setLinearVelocity(objparam.linearVelocity);
       pIt->setRotation(objparam.rotation);
       pIt->setAngularVelocity(objparam.angularVelocity);
-      pIt->setRadius(objparam.radius);
+      pIt->setRadiusAtTemperature(objparam.radiusAtTemperature);
       pIt->setOldContactHistory(objparam.oldContactHistory);
       pIt->setTemperature(objparam.temperature);
 
diff --git a/src/mesa_pd/mpi/notifications/ParticleCopyNotification.h b/src/mesa_pd/mpi/notifications/ParticleCopyNotification.h
index 614f8bea3..b01a4dbeb 100644
--- a/src/mesa_pd/mpi/notifications/ParticleCopyNotification.h
+++ b/src/mesa_pd/mpi/notifications/ParticleCopyNotification.h
@@ -61,7 +61,7 @@ public:
       walberla::mesa_pd::Rot3 rotation {};
       walberla::mesa_pd::Vec3 angularVelocity {real_t(0)};
       walberla::mesa_pd::Vec3 oldTorque {real_t(0)};
-      walberla::real_t radius {real_t(0)};
+      walberla::real_t radiusAtTemperature {real_t(0)};
       uint_t type {0};
       std::map<walberla::id_t, walberla::mesa_pd::data::ContactHistory> oldContactHistory {};
       walberla::real_t temperature {real_t(0)};
@@ -93,7 +93,7 @@ inline data::ParticleStorage::iterator createNewParticle(data::ParticleStorage&
    pIt->setRotation(data.rotation);
    pIt->setAngularVelocity(data.angularVelocity);
    pIt->setOldTorque(data.oldTorque);
-   pIt->setRadius(data.radius);
+   pIt->setRadiusAtTemperature(data.radiusAtTemperature);
    pIt->setType(data.type);
    pIt->setOldContactHistory(data.oldContactHistory);
    pIt->setTemperature(data.temperature);
@@ -140,7 +140,7 @@ mpi::GenericSendBuffer<T,G>& operator<<( mpi::GenericSendBuffer<T,G> & buf, cons
    buf << obj.particle_.getRotation();
    buf << obj.particle_.getAngularVelocity();
    buf << obj.particle_.getOldTorque();
-   buf << obj.particle_.getRadius();
+   buf << obj.particle_.getRadiusAtTemperature();
    buf << obj.particle_.getType();
    buf << obj.particle_.getOldContactHistory();
    buf << obj.particle_.getTemperature();
@@ -168,7 +168,7 @@ mpi::GenericRecvBuffer<T>& operator>>( mpi::GenericRecvBuffer<T> & buf, mesa_pd:
    buf >> objparam.rotation;
    buf >> objparam.angularVelocity;
    buf >> objparam.oldTorque;
-   buf >> objparam.radius;
+   buf >> objparam.radiusAtTemperature;
    buf >> objparam.type;
    buf >> objparam.oldContactHistory;
    buf >> objparam.temperature;
diff --git a/src/mesa_pd/mpi/notifications/ParticleGhostCopyNotification.h b/src/mesa_pd/mpi/notifications/ParticleGhostCopyNotification.h
index aba424550..c4e1c2417 100644
--- a/src/mesa_pd/mpi/notifications/ParticleGhostCopyNotification.h
+++ b/src/mesa_pd/mpi/notifications/ParticleGhostCopyNotification.h
@@ -58,7 +58,7 @@ public:
       size_t shapeID {};
       walberla::mesa_pd::Rot3 rotation {};
       walberla::mesa_pd::Vec3 angularVelocity {real_t(0)};
-      walberla::real_t radius {real_t(0)};
+      walberla::real_t radiusAtTemperature {real_t(0)};
       uint_t type {0};
       std::map<walberla::id_t, walberla::mesa_pd::data::ContactHistory> oldContactHistory {};
       walberla::real_t temperature {real_t(0)};
@@ -83,7 +83,7 @@ inline data::ParticleStorage::iterator createNewParticle(data::ParticleStorage&
    pIt->setShapeID(data.shapeID);
    pIt->setRotation(data.rotation);
    pIt->setAngularVelocity(data.angularVelocity);
-   pIt->setRadius(data.radius);
+   pIt->setRadiusAtTemperature(data.radiusAtTemperature);
    pIt->setType(data.type);
    pIt->setOldContactHistory(data.oldContactHistory);
    pIt->setTemperature(data.temperature);
@@ -123,7 +123,7 @@ mpi::GenericSendBuffer<T,G>& operator<<( mpi::GenericSendBuffer<T,G> & buf, cons
    buf << obj.particle_.getShapeID();
    buf << obj.particle_.getRotation();
    buf << obj.particle_.getAngularVelocity();
-   buf << obj.particle_.getRadius();
+   buf << obj.particle_.getRadiusAtTemperature();
    buf << obj.particle_.getType();
    buf << obj.particle_.getOldContactHistory();
    buf << obj.particle_.getTemperature();
@@ -144,7 +144,7 @@ mpi::GenericRecvBuffer<T>& operator>>( mpi::GenericRecvBuffer<T> & buf, mesa_pd:
    buf >> objparam.shapeID;
    buf >> objparam.rotation;
    buf >> objparam.angularVelocity;
-   buf >> objparam.radius;
+   buf >> objparam.radiusAtTemperature;
    buf >> objparam.type;
    buf >> objparam.oldContactHistory;
    buf >> objparam.temperature;
diff --git a/src/mesa_pd/mpi/notifications/ParticleUpdateNotification.h b/src/mesa_pd/mpi/notifications/ParticleUpdateNotification.h
index 0214aa977..d84508200 100644
--- a/src/mesa_pd/mpi/notifications/ParticleUpdateNotification.h
+++ b/src/mesa_pd/mpi/notifications/ParticleUpdateNotification.h
@@ -50,7 +50,7 @@ public:
    walberla::mesa_pd::Vec3 linearVelocity {real_t(0)};
    walberla::mesa_pd::Rot3 rotation {};
    walberla::mesa_pd::Vec3 angularVelocity {real_t(0)};
-   walberla::real_t radius {real_t(0)};
+   walberla::real_t radiusAtTemperature {real_t(0)};
    std::map<walberla::id_t, walberla::mesa_pd::data::ContactHistory> oldContactHistory {};
    walberla::real_t temperature {real_t(0)};
    };
@@ -87,7 +87,7 @@ mpi::GenericSendBuffer<T,G>& operator<<( mpi::GenericSendBuffer<T,G> & buf, cons
    buf << obj.particle_.getLinearVelocity();
    buf << obj.particle_.getRotation();
    buf << obj.particle_.getAngularVelocity();
-   buf << obj.particle_.getRadius();
+   buf << obj.particle_.getRadiusAtTemperature();
    buf << obj.particle_.getOldContactHistory();
    buf << obj.particle_.getTemperature();
    return buf;
@@ -102,7 +102,7 @@ mpi::GenericRecvBuffer<T>& operator>>( mpi::GenericRecvBuffer<T> & buf, mesa_pd:
    buf >> objparam.linearVelocity;
    buf >> objparam.rotation;
    buf >> objparam.angularVelocity;
-   buf >> objparam.radius;
+   buf >> objparam.radiusAtTemperature;
    buf >> objparam.oldContactHistory;
    buf >> objparam.temperature;
    return buf;
-- 
GitLab