diff --git a/apps/benchmarks/FluidParticleCouplingWithLoadBalancing/Utility.h b/apps/benchmarks/FluidParticleCouplingWithLoadBalancing/Utility.h
index f9e4c05f1273a463486ef7b18647d31020d5bfda..8a394f6bb6ab168cb97841abba5cd54902157845 100644
--- a/apps/benchmarks/FluidParticleCouplingWithLoadBalancing/Utility.h
+++ b/apps/benchmarks/FluidParticleCouplingWithLoadBalancing/Utility.h
@@ -10,21 +10,21 @@ namespace amr {
  * Result from the workload evaluation as described in
  *  Rettinger, Ruede - "Dynamic Load Balancing Techniques for Particulate Flow Simulations", 2019, Computation
  */
-real_t fittedLBMWeightEvaluationFunction(const BlockInfo& blockInfo)
+inline real_t fittedLBMWeightEvaluationFunction(const BlockInfo& blockInfo)
 {
    uint_t Ce = blockInfo.numberOfCells;
    uint_t F  = blockInfo.numberOfFluidCells;
    real_t weight = real_t(7.597476065046571e-06) * real_c(Ce) + real_t(8.95723566283202e-05) * real_c(F) + real_t(-0.1526111388616016);
    return std::max(weight,real_t(0));
 }
-real_t fittedBHWeightEvaluationFunction(const BlockInfo& blockInfo)
+inline real_t fittedBHWeightEvaluationFunction(const BlockInfo& blockInfo)
 {
    uint_t Ce = blockInfo.numberOfCells;
    uint_t NB = blockInfo.numberOfNearBoundaryCells;
    real_t weight = real_t(1.3067711379655123e-07) * real_c(Ce) + real_t(0.0007289549127205142) * real_c(NB) + real_t(-0.1575698071795788);
    return std::max(weight,real_t(0));
 }
-real_t fittedRPDWeightEvaluationFunction(const BlockInfo& blockInfo)
+inline real_t fittedRPDWeightEvaluationFunction(const BlockInfo& blockInfo)
 {
    uint_t Pl = blockInfo.numberOfLocalParticles;
    uint_t Pg = blockInfo.numberOfGhostParticles;
@@ -36,7 +36,7 @@ real_t fittedRPDWeightEvaluationFunction(const BlockInfo& blockInfo)
    real_t weight = real_c(Sc) * ( cPlPg2 * real_c(Pl+Pg) * real_c(Pl+Pg) + cPl * real_c(Pl) + cPg * real_c(Pg) + c );
    return std::max(weight,real_t(0));
 }
-real_t fittedCoup1WeightEvaluationFunction(const BlockInfo& blockInfo)
+inline real_t fittedCoup1WeightEvaluationFunction(const BlockInfo& blockInfo)
 {
    uint_t Ce = blockInfo.numberOfCells;
    uint_t F  = blockInfo.numberOfFluidCells;
@@ -45,7 +45,7 @@ real_t fittedCoup1WeightEvaluationFunction(const BlockInfo& blockInfo)
    real_t weight = real_t(5.610203409278647e-06) * real_c(Ce) + real_t(-7.257635845636656e-07) * real_c(F) + real_t(0.02049703546054693) * real_c(Pl) + real_t(0.04248208493809902) * real_c(Pg) + real_t(-0.26609470510074784);
    return std::max(weight,real_t(0));
 }
-real_t fittedCoup2WeightEvaluationFunction(const BlockInfo& blockInfo)
+inline real_t fittedCoup2WeightEvaluationFunction(const BlockInfo& blockInfo)
 {
    uint_t Ce = blockInfo.numberOfCells;
    uint_t F  = blockInfo.numberOfFluidCells;
@@ -54,7 +54,7 @@ real_t fittedCoup2WeightEvaluationFunction(const BlockInfo& blockInfo)
    real_t weight = real_t(7.198479654682179e-06) * real_c(Ce) + real_t(1.178247475854302e-06) * real_c(F) + real_t(-0.0026401549115124628) * real_c(Pl) + real_t(0.008459646786179298) * real_c(Pg) + real_t(-0.001077320113275954);
    return std::max(weight,real_t(0));
 }
-real_t fittedTotalWeightEvaluationFunction(const BlockInfo& blockInfo)
+inline real_t fittedTotalWeightEvaluationFunction(const BlockInfo& blockInfo)
 {
    return fittedLBMWeightEvaluationFunction(blockInfo) + fittedBHWeightEvaluationFunction(blockInfo) +
           fittedRPDWeightEvaluationFunction(blockInfo) + fittedCoup1WeightEvaluationFunction(blockInfo) +
diff --git a/apps/benchmarks/GranularGas/check.h b/apps/benchmarks/GranularGas/check.h
index 73c457592ce0e78e976e8163e86dc346a16504bc..71652dcdce7fffbf4ee3219150f46f0c8ccfa20c 100644
--- a/apps/benchmarks/GranularGas/check.h
+++ b/apps/benchmarks/GranularGas/check.h
@@ -27,7 +27,7 @@
 namespace walberla {
 namespace mesa_pd {
 
-void check( data::ParticleStorage& ps, blockforest::BlockForest& forest, real_t spacing, const Vec3& shift )
+inline void check( data::ParticleStorage& ps, blockforest::BlockForest& forest, real_t spacing, const Vec3& shift )
 {
    WALBERLA_LOG_INFO_ON_ROOT("*** CHECKING RESULT - START ***");
    auto pIt = ps.begin();
diff --git a/apps/showcases/HeatConduction/ThermalExpansion.h b/apps/showcases/HeatConduction/ThermalExpansion.h
index c56ee0732258220a0b8e67edbaaeb519843d2db3..0818d527d58476684361c9ece4beb5d502611242 100644
--- a/apps/showcases/HeatConduction/ThermalExpansion.h
+++ b/apps/showcases/HeatConduction/ThermalExpansion.h
@@ -47,7 +47,7 @@ private:
    std::vector<real_t> linearExpansionCoefficient_ {};
 };
 
-ThermalExpansion::ThermalExpansion(const uint_t numParticleTypes)
+inline ThermalExpansion::ThermalExpansion(const uint_t numParticleTypes)
 {
    numParticleTypes_ = numParticleTypes;
 
diff --git a/src/lbm/MassEvaluation.h b/src/lbm/MassEvaluation.h
index 7efe2daaebf0d9e89134f77702a889b49e306fea..0a02e1ea065b15cc5cd50a4243e438bc19c9659a 100644
--- a/src/lbm/MassEvaluation.h
+++ b/src/lbm/MassEvaluation.h
@@ -33,7 +33,7 @@ namespace lbm {
 
 namespace internal {
 
-Vector3<real_t> massEvaluationDomain( const shared_ptr< StructuredBlockStorage > & blocks, const uint_t level )
+inline Vector3<real_t> massEvaluationDomain( const shared_ptr< StructuredBlockStorage > & blocks, const uint_t level )
 {
    return Vector3<real_t>( real_c( blocks->getNumberOfXCells(level) ),
                            real_c( blocks->getNumberOfYCells(level) ),
diff --git a/src/lbm/field/initializer/PdfFieldInitializer.impl.h b/src/lbm/field/initializer/PdfFieldInitializer.impl.h
index 757cb180e1e8725b954b75c19a049fbd4139cc86..1a60e350cf171e1c2b435165e23073685e719ed9 100644
--- a/src/lbm/field/initializer/PdfFieldInitializer.impl.h
+++ b/src/lbm/field/initializer/PdfFieldInitializer.impl.h
@@ -25,7 +25,7 @@ namespace lbm {
 namespace initializer {
 
 template<>
-auto getCoordinates<false>(const Cell& globalCell, const real_t dx) {
+inline auto getCoordinates<false>(const Cell& globalCell, const real_t dx) {
 
    Cell coords;
    for(uint_t d = 0; d < 3; ++d) {
@@ -36,7 +36,7 @@ auto getCoordinates<false>(const Cell& globalCell, const real_t dx) {
 }
 
 template<>
-auto getCoordinates<true>(const Cell& globalCell, const real_t dx) {
+inline auto getCoordinates<true>(const Cell& globalCell, const real_t dx) {
 
    Vector3<real_t> coords;
    for(uint_t d = 0; d < 3; ++d) {
diff --git a/src/lbm/free_surface/LoadBalancing.h b/src/lbm/free_surface/LoadBalancing.h
index c9df54478e5c4af223f0be17fe128f86427f21ed..f1638263be4adfd77b7c87d08bfca1444dd52a39 100644
--- a/src/lbm/free_surface/LoadBalancing.h
+++ b/src/lbm/free_surface/LoadBalancing.h
@@ -50,10 +50,10 @@ class ProcessLoadEvaluator;
 /***********************************************************************************************************************
  * Create non-uniform block forest to be used for load balancing.
  **********************************************************************************************************************/
-std::shared_ptr< StructuredBlockForest > createNonUniformBlockForest(const Vector3< uint_t >& domainSize,
-                                                                     const Vector3< uint_t >& cellsPerBlock,
-                                                                     const Vector3< uint_t >& numBlocks,
-                                                                     const Vector3< bool >& periodicity)
+inline std::shared_ptr< StructuredBlockForest > createNonUniformBlockForest(const Vector3< uint_t >& domainSize,
+                                                                            const Vector3< uint_t >& cellsPerBlock,
+                                                                            const Vector3< uint_t >& numBlocks,
+                                                                            const Vector3< bool >& periodicity)
 {
    WALBERLA_CHECK_EQUAL(domainSize[0], cellsPerBlock[0] * numBlocks[0],
                         "The domain size is not divisible by the specified \"cellsPerBlock\" in x-direction.");
diff --git a/src/lbm/free_surface/dynamics/functionality/ReconstructInterfaceCellABB.h b/src/lbm/free_surface/dynamics/functionality/ReconstructInterfaceCellABB.h
index b0f45113c95155be79c6118406bcd911890b9ffa..1172b721d04f9079e58a6c7d68923e8ae3f6f0b3 100644
--- a/src/lbm/free_surface/dynamics/functionality/ReconstructInterfaceCellABB.h
+++ b/src/lbm/free_surface/dynamics/functionality/ReconstructInterfaceCellABB.h
@@ -38,12 +38,12 @@ namespace walberla
 namespace free_surface
 {
 // get index of largest entry in n_dot_ci with isInterfaceOrLiquid==true && isPdfAvailable==false
-uint_t getIndexOfMaximum(const std::vector< bool >& isInterfaceOrLiquid, const std::vector< bool >& isPdfAvailable,
-                         const std::vector< real_t >& n_dot_ci);
+inline uint_t getIndexOfMaximum(const std::vector< bool >& isInterfaceOrLiquid, const std::vector< bool >& isPdfAvailable,
+                                const std::vector< real_t >& n_dot_ci);
 
 // get index of smallest entry in n_dot_ci with isInterfaceOrLiquid==true && isPdfAvailable==false
-uint_t getIndexOfMinimum(const std::vector< bool >& isInterfaceOrLiquid, const std::vector< bool >& isPdfAvailable,
-                         const std::vector< real_t >& n_dot_ci);
+inline uint_t getIndexOfMinimum(const std::vector< bool >& isInterfaceOrLiquid, const std::vector< bool >& isPdfAvailable,
+                                const std::vector< real_t >& n_dot_ci);
 
 // reconstruct PDFs according to pressure anti bounce back boundary condition (page 31, equation 4.5 in dissertation of
 // N. Thuerey, 2007)
@@ -384,8 +384,8 @@ void reconstructInterfaceCellLegacy(const FlagField_T* flagField, const ConstPdf
    }
 }
 
-uint_t getIndexOfMaximum(const std::vector< bool >& isInterfaceOrLiquid, const std::vector< bool >& isPdfAvailable,
-                         const std::vector< real_t >& n_dot_ci)
+inline uint_t getIndexOfMaximum(const std::vector< bool >& isInterfaceOrLiquid, const std::vector< bool >& isPdfAvailable,
+                                const std::vector< real_t >& n_dot_ci)
 {
    real_t maximum = -std::numeric_limits< real_t >::max();
    uint_t index   = std::numeric_limits< uint_t >::max();
@@ -411,8 +411,8 @@ uint_t getIndexOfMaximum(const std::vector< bool >& isInterfaceOrLiquid, const s
    return index;
 }
 
-uint_t getIndexOfMinimum(const std::vector< bool >& isInterfaceOrLiquid, const std::vector< bool >& isPdfAvailable,
-                         const std::vector< real_t >& n_dot_ci)
+inline uint_t getIndexOfMinimum(const std::vector< bool >& isInterfaceOrLiquid, const std::vector< bool >& isPdfAvailable,
+                                const std::vector< real_t >& n_dot_ci)
 {
    real_t minimum = std::numeric_limits< real_t >::max();
    uint_t index   = std::numeric_limits< uint_t >::max();
diff --git a/src/lbm_mesapd_coupling/amr/InfoCollection.h b/src/lbm_mesapd_coupling/amr/InfoCollection.h
index 621ab8f7b2826d8d3a4ce3ebbb5d9818cfc87fe9..8783a945cdd9c3831e2b295c5f977209f2de6111 100644
--- a/src/lbm_mesapd_coupling/amr/InfoCollection.h
+++ b/src/lbm_mesapd_coupling/amr/InfoCollection.h
@@ -115,6 +115,7 @@ void updateAndSyncInfoCollection(BlockForest& bf, const BlockDataID boundaryHand
    }
 }
 
+inline
 void getBlockInfoFromInfoCollection( const PhantomBlock * block, const shared_ptr<InfoCollection>& ic,
                                      BlockInfo & blockInfo )
 {
diff --git a/src/lbm_mesapd_coupling/mapping/ParticleBoundingBox.h b/src/lbm_mesapd_coupling/mapping/ParticleBoundingBox.h
index aa32158a99094b283016cace0a35632aae300dff..7c2f49a3b8ea4808d9ea8306917579c83b0c3d52 100644
--- a/src/lbm_mesapd_coupling/mapping/ParticleBoundingBox.h
+++ b/src/lbm_mesapd_coupling/mapping/ParticleBoundingBox.h
@@ -36,9 +36,9 @@ namespace lbm_mesapd_coupling {
  * Obtain a block-local cell bounding box from a given AABB (e.g. the particle's AABB)
  * If the given AABB is (partly) infinite, AABBIsInfinite should be set to true (e.g. for infinite particles)
  */
-CellInterval getCellBBFromAABB( const math::AABB & aabb, bool AABBIsInfinite,
-                                const IBlock & block, StructuredBlockStorage & blockStorage,
-                                const uint_t numberOfGhostLayersToInclude)
+inline CellInterval getCellBBFromAABB( const math::AABB & aabb, bool AABBIsInfinite,
+                                       const IBlock & block, StructuredBlockStorage & blockStorage,
+                                       const uint_t numberOfGhostLayersToInclude)
 {
 
    CellInterval cellBB;
diff --git a/src/lbm_mesapd_coupling/utility/OmegaBulkAdaption.h b/src/lbm_mesapd_coupling/utility/OmegaBulkAdaption.h
index 226c03c2f16a2a8ba1b399545a5f38057ed29340..f8a60c3a943e3e046d45038c8d86513f5c088c55 100644
--- a/src/lbm_mesapd_coupling/utility/OmegaBulkAdaption.h
+++ b/src/lbm_mesapd_coupling/utility/OmegaBulkAdaption.h
@@ -29,20 +29,20 @@ namespace lbm_mesapd_coupling {
 
 // utility functions
 
-real_t bulkViscosityFromOmegaBulk(real_t omegaBulk)
+inline real_t bulkViscosityFromOmegaBulk(real_t omegaBulk)
 {
    return real_t(2) / real_t(9) * ( real_t(1) / omegaBulk - real_t(0.5) );
 }
 
 
-real_t omegaBulkFromBulkViscosity(real_t bulkViscosity)
+inline real_t omegaBulkFromBulkViscosity(real_t bulkViscosity)
 {
    return real_t(2) / ( real_t(9) * bulkViscosity + real_t(1) );
 }
 
 // see Khirevich et al. - Coarse- and fine-grid numerical behavior of MRT/TRT lattice-Boltzmann schemes in regular and random sphere packings
 // LambdaBulk is the "magic parameter" here, i.e. the ratio between Lambda_e and Lambda_nu, Eq. 19
-real_t omegaBulkFromOmega(real_t omega, real_t LambdaBulk = real_t(1))
+inline real_t omegaBulkFromOmega(real_t omega, real_t LambdaBulk = real_t(1))
 {
    return real_t(1) / (LambdaBulk * ( real_t(1) / omega - real_t(1)/ real_t(2) ) + real_t(1)/ real_t(2) );
 }
diff --git a/src/mesa_pd/collision_detection/GeneralContactDetection.h b/src/mesa_pd/collision_detection/GeneralContactDetection.h
index 9564916c02a97432793a09f5805bf9f9072f73ac..c9101312667fde8afaa7761a05201dde518fec97 100644
--- a/src/mesa_pd/collision_detection/GeneralContactDetection.h
+++ b/src/mesa_pd/collision_detection/GeneralContactDetection.h
@@ -191,7 +191,7 @@ bool GeneralContactDetection::operator()(const size_t idx1,
    return operator()(idx2, idx1, geo2, geo1, ac);
 }
 
-bool GeneralContactDetection::collideGJKEPA(Support& geom0, Support& geom1)
+inline bool GeneralContactDetection::collideGJKEPA(Support& geom0, Support& geom1)
 {
    real_t margin = real_t(1e-4);
    GJK gjk;
diff --git a/src/mesa_pd/common/AABBConversion.h b/src/mesa_pd/common/AABBConversion.h
index c105d1ed2b336b91c671e203bba786601b613ffd..be579b9701d74a4879046327f331c0c5131d3b7d 100644
--- a/src/mesa_pd/common/AABBConversion.h
+++ b/src/mesa_pd/common/AABBConversion.h
@@ -30,7 +30,7 @@
 namespace walberla {
 namespace mesa_pd {
 
-math::AABB getAABBFromInteractionRadius(const Vector3<real_t> & pos, const real_t interactionRadius )
+inline math::AABB getAABBFromInteractionRadius(const Vector3<real_t> & pos, const real_t interactionRadius )
 {
    WALBERLA_ASSERT_GREATER(interactionRadius, 0_r, "Did you forget to set the interaction radius?");
    return math::AABB( pos[0]-interactionRadius, pos[1]-interactionRadius, pos[2]-interactionRadius,
diff --git a/src/mesa_pd/common/Contains.h b/src/mesa_pd/common/Contains.h
index 00e74bb3990d29436d1a4507796d315c5d050ea7..8c707388bda2b541dee4eecab28518fac685d8a5 100644
--- a/src/mesa_pd/common/Contains.h
+++ b/src/mesa_pd/common/Contains.h
@@ -42,18 +42,21 @@ namespace mesa_pd {
  * or in body frame coordinates (BF) which requires the point to be first transformed
  */
 
+inline
 bool isPointInsideSphere(const Vec3& point,
                          const Vec3& spherePosition, const real_t sphereRadius )
 {
    return !((point - spherePosition).sqrLength() > sphereRadius * sphereRadius);
 }
 
+inline
 bool isPointInsideHalfSpace(const Vec3& point,
                             const Vec3& halfSpacePosition, const Vec3& halfSpaceNormal )
 {
    return !((point - halfSpacePosition) * halfSpaceNormal > real_t(0));
 }
 
+inline
 bool isPointInsideBoxBF(const Vec3& pointBF,
                         const Vec3& edgeLengths )
 {
@@ -62,6 +65,7 @@ bool isPointInsideBoxBF(const Vec3& pointBF,
           std::fabs(pointBF[2]) <= real_t(0.5)*edgeLengths[2];
 }
 
+inline
 bool isPointInsideEllipsoidBF(const Vec3& pointBF,
                               const Vec3& semiAxes )
 {
@@ -69,6 +73,7 @@ bool isPointInsideEllipsoidBF(const Vec3& pointBF,
             + (pointBF[2] * pointBF[2])/(semiAxes[2] * semiAxes[2]) <= 1_r );
 }
 
+inline
 bool isPointInsideCylindricalBoundary(const Vec3& point,
                                       const Vec3& cylindricalBoundaryPosition, const real_t radius, const Vec3& axis  )
 {
@@ -77,6 +82,7 @@ bool isPointInsideCylindricalBoundary(const Vec3& point,
 }
 
 #ifdef WALBERLA_MESAPD_CONVEX_POLYHEDRON_AVAILABLE
+inline
 bool isPointInsideConvexPolyhedronBF(const Vec3& point, const mesh::TriangleMesh& mesh)
 {
    WALBERLA_ASSERT(mesh.has_face_normals(), "Provided mesh has no face normals! E.g., call `mesh.request_face_normals(); mesh.update_face_normals();` to add them.")
diff --git a/src/mesa_pd/common/RayParticleIntersection.h b/src/mesa_pd/common/RayParticleIntersection.h
index b62306aa5dd6014d301ddbd6943a8a0f45b700b9..79c7910a6a3aae6108383a8b8e7b0729dd228d20 100644
--- a/src/mesa_pd/common/RayParticleIntersection.h
+++ b/src/mesa_pd/common/RayParticleIntersection.h
@@ -38,6 +38,7 @@ namespace mesa_pd {
  * or in body frame coordinates (BF) which requires the point to be first transformed
  */
 
+inline
 real_t raySphereIntersectionRatio( const Vec3& rayOrigin, const Vec3& rayDirection,
                                    const Vec3& spherePosition, const real_t sphereRadius )
 {
@@ -62,6 +63,7 @@ real_t raySphereIntersectionRatio( const Vec3& rayOrigin, const Vec3& rayDirecti
    return delta;
 }
 
+inline
 real_t rayHalfSpaceIntersectionRatio( const Vec3& rayOrigin, const Vec3& rayDirection,
                                       const Vec3& halfSpacePosition, const Vec3& halfSpaceNormal)
 {
@@ -82,6 +84,7 @@ real_t rayHalfSpaceIntersectionRatio( const Vec3& rayOrigin, const Vec3& rayDire
    return delta;
 }
 
+inline
 real_t rayEllipsoidIntersectionRatioBF( const Vec3& rayOriginBF, const Vec3& rayDirectionBF,
                                         const Vec3& ellipsoidSemiAxes)
 {
diff --git a/src/mesa_pd/data/HashGrids.h b/src/mesa_pd/data/HashGrids.h
index 83c210b6534a5a1b0fba9675a1585d29f8522486..d847aef68b72f1aec5c0a9d6584b5f3945fa0000 100644
--- a/src/mesa_pd/data/HashGrids.h
+++ b/src/mesa_pd/data/HashGrids.h
@@ -351,7 +351,7 @@ private:
 
 // HashGrids::HashGrid member function implementations
 
-HashGrids::HashGrid::HashGrid( real_t cellSpan )
+inline HashGrids::HashGrid::HashGrid( real_t cellSpan )
 {
    // Initialization of all member variables and ...
    xCellCount_   = math::uintIsPowerOfTwo( xCellCount ) ? xCellCount : 16;
@@ -386,7 +386,7 @@ HashGrids::HashGrid::HashGrid( real_t cellSpan )
    particleCount_ = 0;
 }
 
-HashGrids::HashGrid::~HashGrid()
+inline HashGrids::HashGrid::~HashGrid()
 {
    clear();
 
@@ -396,7 +396,7 @@ HashGrids::HashGrid::~HashGrid()
    delete[] cell_;
 }
 
-void HashGrids::HashGrid::clear()
+inline void HashGrids::HashGrid::clear()
 {
    for( auto cellIt = occupiedCells_.begin(); cellIt < occupiedCells_.end(); ++cellIt ) {
       delete (*cellIt)->particles_;
@@ -412,7 +412,7 @@ void HashGrids::HashGrid::clear()
  * This function is used to initialize the offset arrays of all grid cells. The offsets are required
  * for ensuring fast direct access to all directly adjacent cells for each cell in the hash grid.
  */
-void HashGrids::HashGrid::initializeNeighborOffsets()
+inline void HashGrids::HashGrid::initializeNeighborOffsets()
 {
    offset_t xc   = static_cast<offset_t>( xCellCount_ );
    offset_t yc   = static_cast<offset_t>( yCellCount_ );
@@ -516,7 +516,7 @@ void HashGrids::HashGrid::addParticle( size_t p_idx, Accessor& ac )
 //*************************************************************************************************
 /*!\brief Adds a particle to a specific cell in this hash grid.
  */
-void HashGrids::HashGrid::addParticleToCell( size_t p_idx, Cell* cell )
+inline void HashGrids::HashGrid::addParticleToCell( size_t p_idx, Cell* cell )
 {
    // If this cell is already occupied by other particles, which means the pointer to the particle
    // container holds a valid address and thus the container itself is properly initialized, then
@@ -564,7 +564,7 @@ size_t HashGrids::HashGrid::hashOfParticle( size_t p_idx, Accessor& ac ) const
  * Note that the modulo calculations are replaced with fast bitwise AND operations - hence, the
  * spatial dimensions of the hash grid must be restricted to powers of two!
  */
-size_t HashGrids::HashGrid::hashPoint(real_t x, real_t y, real_t z) const {
+inline size_t HashGrids::HashGrid::hashPoint(real_t x, real_t y, real_t z) const {
    size_t xHash;
    size_t yHash;
    size_t zHash;
@@ -773,7 +773,7 @@ void HashGrids::HashGrid::checkAgainstVectorEachParticlePairHalf( const Particle
 //*************************************************************************************************
 
 // clear all particles and grids
-void HashGrids::clearAll()
+inline void HashGrids::clearAll()
 {
    gridList_.clear();
    infiniteParticles_.clear();
@@ -781,7 +781,7 @@ void HashGrids::clearAll()
 
 // clear only all particles from the hash grids, but maintain the overall grid hierarchy
 // useful for "updating" the data structure in each time step, but also prevents clean-up of unnecessary grids
-void HashGrids::clear()
+inline void HashGrids::clear()
 {
    for( auto gridIt = gridList_.begin(); gridIt != gridList_.end(); ++gridIt ) {
       (*gridIt)->clear();
@@ -849,7 +849,7 @@ void HashGrids::operator()(const size_t p_idx, Accessor& ac)
    }
 }
 
-void HashGrids::addInfiniteParticle(size_t p_idx)
+inline void HashGrids::addInfiniteParticle(size_t p_idx)
 {
    infiniteParticles_.push_back(p_idx);
 }
diff --git a/src/mesa_pd/data/LinkedCells.h b/src/mesa_pd/data/LinkedCells.h
index 9ba699fd0da4a36afc1d76472b1c37c1917a3bde..26ac5c7972edca90e5c5e4833ae42d7d8ff152a6 100644
--- a/src/mesa_pd/data/LinkedCells.h
+++ b/src/mesa_pd/data/LinkedCells.h
@@ -171,7 +171,7 @@ LinkedCells::LinkedCells(const math::AABB& domain, const Vec3& cellDiameter)
    std::fill(cells_.begin(), cells_.end(), -1);
 }
 
-void LinkedCells::clear()
+inline void LinkedCells::clear()
 {
    const uint64_t cellsSize = cells_.size();
    //clear existing linked cells
diff --git a/src/mesa_pd/data/SparseLinkedCells.h b/src/mesa_pd/data/SparseLinkedCells.h
index 2c184de476a67c8a02eaf19aae72540d625e2bd9..3e98b174c3a8c2a0c2985891259e0463f902c198 100644
--- a/src/mesa_pd/data/SparseLinkedCells.h
+++ b/src/mesa_pd/data/SparseLinkedCells.h
@@ -187,7 +187,7 @@ SparseLinkedCells::SparseLinkedCells(const math::AABB& domain, const Vec3& cellD
    std::fill(cells_.begin(), cells_.end(), -1);
 }
 
-void SparseLinkedCells::clear()
+inline void SparseLinkedCells::clear()
 {
    for (const auto v : nonEmptyCells_)
    {
diff --git a/src/mesa_pd/kernel/ForceLJ.h b/src/mesa_pd/kernel/ForceLJ.h
index 1affbc8f82ca99fa1fae40e2fb3fef1f22cbdd58..14d0877728fc578e67a7021639c2e404827cd682 100644
--- a/src/mesa_pd/kernel/ForceLJ.h
+++ b/src/mesa_pd/kernel/ForceLJ.h
@@ -80,7 +80,7 @@ private:
    std::vector<real_t> sigma {};
 };
 
-ForceLJ::ForceLJ(const uint_t numParticleTypes)
+inline ForceLJ::ForceLJ(const uint_t numParticleTypes)
 {
    numParticleTypes_ = numParticleTypes;
    
diff --git a/src/mesa_pd/kernel/HeatConduction.h b/src/mesa_pd/kernel/HeatConduction.h
index e41d117047bce400a4e9d6601c70aa1d46fe7f69..efedff4efaf8b6bba407a141f13155ae4e294c9d 100644
--- a/src/mesa_pd/kernel/HeatConduction.h
+++ b/src/mesa_pd/kernel/HeatConduction.h
@@ -81,7 +81,7 @@ private:
    std::vector<real_t> conductance_ {};
 };
 
-HeatConduction::HeatConduction(const uint_t numParticleTypes)
+inline HeatConduction::HeatConduction(const uint_t numParticleTypes)
 {
    numParticleTypes_ = numParticleTypes;
    
diff --git a/src/mesa_pd/kernel/LinearSpringDashpot.h b/src/mesa_pd/kernel/LinearSpringDashpot.h
index b4b35082cba39c5884b757a722baa494a5ad6750..db9ca2e902bbf53ebba37a850875a5cf04bde4a8 100644
--- a/src/mesa_pd/kernel/LinearSpringDashpot.h
+++ b/src/mesa_pd/kernel/LinearSpringDashpot.h
@@ -133,7 +133,7 @@ private:
    std::vector<real_t> frictionCoefficientDynamic_ {};
 };
 
-LinearSpringDashpot::LinearSpringDashpot(const uint_t numParticleTypes)
+inline LinearSpringDashpot::LinearSpringDashpot(const uint_t numParticleTypes)
 {
    numParticleTypes_ = numParticleTypes;
    
diff --git a/src/mesa_pd/kernel/NonLinearSpringDashpot.h b/src/mesa_pd/kernel/NonLinearSpringDashpot.h
index b5d732c6d2b27ea535c62404baa3b75d81ace8a6..fc1e0515be1d3f7c90dbb6f30744ab1f729a75e5 100644
--- a/src/mesa_pd/kernel/NonLinearSpringDashpot.h
+++ b/src/mesa_pd/kernel/NonLinearSpringDashpot.h
@@ -120,7 +120,7 @@ private:
    std::vector<real_t> frictionCoefficientDynamic_ {};
 };
 
-NonLinearSpringDashpot::NonLinearSpringDashpot(const uint_t numParticleTypes, const real_t collisionTime)
+inline NonLinearSpringDashpot::NonLinearSpringDashpot(const uint_t numParticleTypes, const real_t collisionTime)
 {
    numParticleTypes_ = numParticleTypes;
    
diff --git a/src/mesa_pd/kernel/SpringDashpot.h b/src/mesa_pd/kernel/SpringDashpot.h
index 7c59af5dcd977c56ff937a5f0aeb46e5f7f557ba..eecd2ee0f302db8efddad1986b35a768b40db391 100644
--- a/src/mesa_pd/kernel/SpringDashpot.h
+++ b/src/mesa_pd/kernel/SpringDashpot.h
@@ -132,7 +132,7 @@ private:
    std::vector<real_t> friction_ {};
 };
 
-SpringDashpot::SpringDashpot(const uint_t numParticleTypes)
+inline SpringDashpot::SpringDashpot(const uint_t numParticleTypes)
 {
    numParticleTypes_ = numParticleTypes;
    
diff --git a/src/mesa_pd/kernel/SpringDashpotSpring.h b/src/mesa_pd/kernel/SpringDashpotSpring.h
index 07eed81524cb9712b3ee027886947b28c973e0ce..ece97ace1ab459e82a21050e7a356df402c424d0 100644
--- a/src/mesa_pd/kernel/SpringDashpotSpring.h
+++ b/src/mesa_pd/kernel/SpringDashpotSpring.h
@@ -116,7 +116,7 @@ private:
    std::vector<real_t> coefficientOfFriction_ {};
 };
 
-SpringDashpotSpring::SpringDashpotSpring(const uint_t numParticleTypes)
+inline SpringDashpotSpring::SpringDashpotSpring(const uint_t numParticleTypes)
 {
    numParticleTypes_ = numParticleTypes;
    
diff --git a/src/mesa_pd/kernel/TemperatureIntegration.h b/src/mesa_pd/kernel/TemperatureIntegration.h
index 8f8a771cec0217ab8b4d096637e50c850cb8ba36..905973f8da70828b7e68de953c759ad8df4f21cc 100644
--- a/src/mesa_pd/kernel/TemperatureIntegration.h
+++ b/src/mesa_pd/kernel/TemperatureIntegration.h
@@ -81,7 +81,7 @@ private:
    std::vector<real_t> invSpecificHeat_ {};
 };
 
-TemperatureIntegration::TemperatureIntegration(const real_t dt, const uint_t numParticleTypes)
+inline TemperatureIntegration::TemperatureIntegration(const real_t dt, const uint_t numParticleTypes)
    : dt_(dt)
 {
    numParticleTypes_ = numParticleTypes;
diff --git a/src/mesa_pd/kernel/cnt/IsotropicVDWContact.h b/src/mesa_pd/kernel/cnt/IsotropicVDWContact.h
index a92ecda09637c90da12722611665836875016328..c3ffbfffe80be84e453bd5c51356a2d222718b31 100644
--- a/src/mesa_pd/kernel/cnt/IsotropicVDWContact.h
+++ b/src/mesa_pd/kernel/cnt/IsotropicVDWContact.h
@@ -95,7 +95,7 @@ void IsotropicVDWContact::operator()(const size_t p_idx1,
    addForceAtomic(p_idx2, ac, -force);
 }
 
-real_t IsotropicVDWContact::equilibriumDistance()
+inline real_t IsotropicVDWContact::equilibriumDistance()
 {
    return r * ( std::pow( (alpha*A)/(beta*B), 1_r/(alpha-beta)) + 2_r);
 }
diff --git a/src/mesa_pd/mpi/ReduceContactHistory.h b/src/mesa_pd/mpi/ReduceContactHistory.h
index efd4e9bd85dfe65d2d38c9be5bb9fac4250ed652..a741c79526874f30fda0d5b90cc5b4555ed77e49 100644
--- a/src/mesa_pd/mpi/ReduceContactHistory.h
+++ b/src/mesa_pd/mpi/ReduceContactHistory.h
@@ -85,7 +85,7 @@ private:
    int numProcesses_ = walberla::mpi::MPIManager::instance()->numProcesses();
 };
 
-void ReduceContactHistory::operator()(data::ParticleStorage& ps) const
+inline void ReduceContactHistory::operator()(data::ParticleStorage& ps) const
 {
    //no need to reduce if run with only one process
    if (numProcesses_ != 1)
diff --git a/src/mesa_pd/mpi/notifications/ContactHistoryNotification.h b/src/mesa_pd/mpi/notifications/ContactHistoryNotification.h
index 3785f0ffbbb0e92e184c9635d1626e95ec9ef388..7d46c26de144724624bb164412695335436d4f82 100644
--- a/src/mesa_pd/mpi/notifications/ContactHistoryNotification.h
+++ b/src/mesa_pd/mpi/notifications/ContactHistoryNotification.h
@@ -58,12 +58,12 @@ public:
 };
 
 template <>
-void reset<ContactHistoryNotification>(data::Particle& p)
+inline void reset<ContactHistoryNotification>(data::Particle& p)
 {
    p.setNewContactHistory(std::map<walberla::id_t, walberla::mesa_pd::data::ContactHistory>());
 }
 
-void reduce(data::Particle&& p, const ContactHistoryNotification::Parameters& objparam)
+inline void reduce(data::Particle&& p, const ContactHistoryNotification::Parameters& objparam)
 {
    auto& ch = p.getNewContactHistoryRef();
    for (auto& entry : objparam.contactHistory_)
diff --git a/src/mesa_pd/mpi/notifications/ForceTorqueNotification.h b/src/mesa_pd/mpi/notifications/ForceTorqueNotification.h
index 5c8ad6f79507d114423115d9e3082808a166e3e1..54fe3902675020ab1fd9d5f2a61e3fe3b40bbc03 100644
--- a/src/mesa_pd/mpi/notifications/ForceTorqueNotification.h
+++ b/src/mesa_pd/mpi/notifications/ForceTorqueNotification.h
@@ -57,19 +57,19 @@ public:
 };
 
 template <>
-void reset<ForceTorqueNotification>(data::Particle& p)
+inline void reset<ForceTorqueNotification>(data::Particle& p)
 {
    p.setForce( Vec3(real_t(0)) );
    p.setTorque( Vec3(real_t(0)) );
 }
 
-void reduce(data::Particle&& p, const ForceTorqueNotification::Parameters& objparam)
+inline void reduce(data::Particle&& p, const ForceTorqueNotification::Parameters& objparam)
 {
    p.getForceRef() += objparam.force_;
    p.getTorqueRef() += objparam.torque_;
 }
 
-void update(data::Particle&& p, const ForceTorqueNotification::Parameters& objparam)
+inline void update(data::Particle&& p, const ForceTorqueNotification::Parameters& objparam)
 {
    p.setForce( objparam.force_ );
    p.setTorque( objparam.torque_ );
diff --git a/src/mesa_pd/mpi/notifications/HeatFluxNotification.h b/src/mesa_pd/mpi/notifications/HeatFluxNotification.h
index 6b580a50083295eb7c40646ef95dbc5af94c1d34..53ebde522e2a70606e0af6e1471914ecac692426 100644
--- a/src/mesa_pd/mpi/notifications/HeatFluxNotification.h
+++ b/src/mesa_pd/mpi/notifications/HeatFluxNotification.h
@@ -56,17 +56,17 @@ public:
 };
 
 template <>
-void reset<HeatFluxNotification>(data::Particle& p)
+inline void reset<HeatFluxNotification>(data::Particle& p)
 {
    p.setHeatFlux( real_t(0) );
 }
 
-void reduce(data::Particle&& p, const HeatFluxNotification::Parameters& objparam)
+inline void reduce(data::Particle&& p, const HeatFluxNotification::Parameters& objparam)
 {
    p.getHeatFluxRef() += objparam.heatFlux_;
 }
 
-void update(data::Particle&& p, const HeatFluxNotification::Parameters& objparam)
+inline void update(data::Particle&& p, const HeatFluxNotification::Parameters& objparam)
 {
    p.setHeatFlux( objparam.heatFlux_ );
 }
diff --git a/src/mesa_pd/mpi/notifications/HydrodynamicForceTorqueNotification.h b/src/mesa_pd/mpi/notifications/HydrodynamicForceTorqueNotification.h
index ce874188127bd4cfa0c11e6d7087ce7e1ef61bd5..a8f1c3376f9a5af5ce4cee6b45d075ad78846553 100644
--- a/src/mesa_pd/mpi/notifications/HydrodynamicForceTorqueNotification.h
+++ b/src/mesa_pd/mpi/notifications/HydrodynamicForceTorqueNotification.h
@@ -57,19 +57,19 @@ public:
 };
 
 template <>
-void reset<HydrodynamicForceTorqueNotification>(data::Particle& p)
+inline void reset<HydrodynamicForceTorqueNotification>(data::Particle& p)
 {
    p.setHydrodynamicForce( Vec3(real_t(0)) );
    p.setHydrodynamicTorque( Vec3(real_t(0)) );
 }
 
-void reduce(data::Particle&& p, const HydrodynamicForceTorqueNotification::Parameters& objparam)
+inline void reduce(data::Particle&& p, const HydrodynamicForceTorqueNotification::Parameters& objparam)
 {
    p.getHydrodynamicForceRef() += objparam.hydrodynamicForce_;
    p.getHydrodynamicTorqueRef() += objparam.hydrodynamicTorque_;
 }
 
-void update(data::Particle&& p, const HydrodynamicForceTorqueNotification::Parameters& objparam)
+inline void update(data::Particle&& p, const HydrodynamicForceTorqueNotification::Parameters& objparam)
 {
    p.setHydrodynamicForce( objparam.hydrodynamicForce_ );
    p.setHydrodynamicTorque( objparam.hydrodynamicTorque_ );
diff --git a/src/mesa_pd/mpi/notifications/VelocityCorrectionNotification.h b/src/mesa_pd/mpi/notifications/VelocityCorrectionNotification.h
index 7f9347f918ca77ea079c12988e8d3d03f3254fea..6d4a88399d6655a361b334f5599d4effae83df27 100644
--- a/src/mesa_pd/mpi/notifications/VelocityCorrectionNotification.h
+++ b/src/mesa_pd/mpi/notifications/VelocityCorrectionNotification.h
@@ -61,14 +61,14 @@ public:
 
 
 // Reduce method for reduction (add up the velocity corrections)
-void reduce(data::Particle&& p, const VelocityCorrectionNotification::Parameters& objparam)
+inline void reduce(data::Particle&& p, const VelocityCorrectionNotification::Parameters& objparam)
 {
    p.getDvRef() += objparam.dv_;
    p.getDwRef() += objparam.dw_;
 }
 
 template<>
-void reset<VelocityCorrectionNotification>(data::Particle& p )
+inline void reset<VelocityCorrectionNotification>(data::Particle& p )
 {
    p.setDv( Vec3(real_t(0)) );
    p.setDw( Vec3(real_t(0)) );
diff --git a/src/mesa_pd/mpi/notifications/VelocityUpdateNotification.h b/src/mesa_pd/mpi/notifications/VelocityUpdateNotification.h
index f526d99c8f21c50f7b2b7ac224498d3f58d19125..7c2ef66d63079342f4e3c050cf39d1b4c5a18a8a 100644
--- a/src/mesa_pd/mpi/notifications/VelocityUpdateNotification.h
+++ b/src/mesa_pd/mpi/notifications/VelocityUpdateNotification.h
@@ -66,7 +66,7 @@ public:
 real_t VelocityUpdateNotification::Parameters::relaxationParam = real_t(0.8);
 
 // Update method for broadcast
-void update(data::Particle&& p, const VelocityUpdateNotification::Parameters& objparam) {
+inline void update(data::Particle&& p, const VelocityUpdateNotification::Parameters& objparam) {
    // Reset the velocity corrections dv/dw of ghost particle
    p.getDvRef() = Vec3();
    p.getDwRef() = Vec3();
@@ -75,7 +75,7 @@ void update(data::Particle&& p, const VelocityUpdateNotification::Parameters& ob
 }
 
 template<>
-void reset<VelocityUpdateNotification>(data::Particle& p )
+inline void reset<VelocityUpdateNotification>(data::Particle& p )
 {
    p.setDv( Vec3(real_t(0)) );
    p.setDw( Vec3(real_t(0)) );
diff --git a/src/pde/sweeps/Multigrid.impl.h b/src/pde/sweeps/Multigrid.impl.h
index da7c0ab6615781b492adcb4e6a9ce55c804e0d6f..47080ea6385542720405ff4ab257dc2331bce38e 100644
--- a/src/pde/sweeps/Multigrid.impl.h
+++ b/src/pde/sweeps/Multigrid.impl.h
@@ -168,7 +168,7 @@ void CoarsenStencilFieldsDCA<Stencil_T >::operator()( const std::vector<BlockDat
 
 
 template< >
-void CoarsenStencilFieldsGCA< stencil::D3Q7 >::operator()( const std::vector<BlockDataID> & stencilFieldId ) const
+inline void CoarsenStencilFieldsGCA< stencil::D3Q7 >::operator()( const std::vector<BlockDataID> & stencilFieldId ) const
 {
 
    WALBERLA_ASSERT_EQUAL(numLvl_, stencilFieldId.size(), "This function can only be called when operating with stencil fields!");
diff --git a/src/pe_coupling/discrete_particle_methods/correlations/AddedMassForceCorrelations.h b/src/pe_coupling/discrete_particle_methods/correlations/AddedMassForceCorrelations.h
index 4574cbc19d02d50b964cb402ba2e8595f5fb36f1..81651697744351f57a7597600f74211a3dd74cbe 100644
--- a/src/pe_coupling/discrete_particle_methods/correlations/AddedMassForceCorrelations.h
+++ b/src/pe_coupling/discrete_particle_methods/correlations/AddedMassForceCorrelations.h
@@ -36,6 +36,7 @@ namespace discrete_particle_methods {
  *                   const real_t & bodyVolume, const real_t & fluidDensity )
  */
 
+inline
 Vector3<real_t> addedMassForceFinn( const Vector3<real_t> & timeDerivativeFluidVel, const Vector3<real_t> & timeDerivativeBodyVel,
                                     const real_t & bodyVolume, const real_t & fluidDensity )
 {
@@ -44,6 +45,7 @@ Vector3<real_t> addedMassForceFinn( const Vector3<real_t> & timeDerivativeFluidV
    return bodyVolume * fluidDensity * Coeffam * ( timeDerivativeFluidVel - timeDerivativeBodyVel );
 }
 
+inline
 Vector3<real_t> noAddedMassForce( const Vector3<real_t> &, const Vector3<real_t> &, const real_t &, const real_t & )
 {
    return Vector3<real_t>(real_t(0));
diff --git a/src/pe_coupling/discrete_particle_methods/correlations/DragForceCorrelations.h b/src/pe_coupling/discrete_particle_methods/correlations/DragForceCorrelations.h
index 1db90c9457a398846ba1f478ee33ebc3a5da2333..002be4a30685445cd5806423dec90b721c134be1 100644
--- a/src/pe_coupling/discrete_particle_methods/correlations/DragForceCorrelations.h
+++ b/src/pe_coupling/discrete_particle_methods/correlations/DragForceCorrelations.h
@@ -45,7 +45,7 @@ namespace discrete_particle_methods {
 
 // equation to calculate the drag coefficient on isolated spherical particle
 // Schiller, L., Naumann, A., 1935. A drag coefficient correlation. Vdi Zeitung 77, 318-320.
-real_t dragCoeffSchillerNaumann( real_t reynoldsNumber )
+inline real_t dragCoeffSchillerNaumann( real_t reynoldsNumber )
 {
    WALBERLA_ASSERT_GREATER_EQUAL( reynoldsNumber, real_t(0) );
 
@@ -55,7 +55,7 @@ real_t dragCoeffSchillerNaumann( real_t reynoldsNumber )
 
 // Coefficient from Stokes' law for drag, only valid for Stokes regime (low Reynolds numbers)
 // = 3 * math::pi * mu * D * fluidVolumeFraction
-real_t dragCoeffStokes ( real_t fluidVolumeFraction, real_t diameter, real_t fluidDynamicViscosity )
+inline real_t dragCoeffStokes ( real_t fluidVolumeFraction, real_t diameter, real_t fluidDynamicViscosity )
 {
    return real_t(3) * math::pi * diameter * fluidDynamicViscosity * fluidVolumeFraction;
 }
@@ -72,6 +72,7 @@ const real_t thresholdAbsoluteVelocityDifference = real_t(1e-10);
 //////////////////////
 
 // Stokes drag law
+inline
 Vector3<real_t> dragForceStokes( const Vector3<real_t> & fluidVel, const Vector3<real_t> & particleVel,
                                  real_t solidVolumeFraction, real_t diameter, real_t fluidDynamicViscosity, real_t /*fluidDensity*/ )
 {
@@ -92,6 +93,7 @@ Vector3<real_t> dragForceStokes( const Vector3<real_t> & fluidVel, const Vector3
 // S. Ergun, Fluid flow through packed columns. Chemical Engineering Progress 48 (1952), 89-94.
 // Y. C. Wen, Y.H. Yu, Mechanics of fluidization. Chemical Engineering Progress Symposium Series 62 (1966), 100-111.
 // see also Beetstra, van der Hoef, Kuipers, "Drag Force of Intermediate Reynolds Number Flow Past Mono- and Bidisperse Arrays of Spheres" (2007)
+inline
 Vector3<real_t>  dragForceErgunWenYu( const Vector3<real_t> & fluidVel, const Vector3<real_t> & particleVel,
                                       real_t solidVolumeFraction, real_t diameter, real_t fluidDynamicViscosity, real_t fluidDensity )
 {
@@ -123,6 +125,7 @@ Vector3<real_t>  dragForceErgunWenYu( const Vector3<real_t> & fluidVel, const Ve
 
 // drag correlation proposed by Tang et al. - "A New Drag Correlation from Fully Resolved Simulations of Flow Past
 // Monodisperse Static Arrays of Spheres", AiChE, 2014
+inline
 Vector3<real_t> dragForceTang( const Vector3<real_t> & fluidVel, const Vector3<real_t> & particleVel,
                                real_t solidVolumeFraction, real_t diameter, real_t fluidDynamicViscosity, real_t fluidDensity )
 {
@@ -151,6 +154,7 @@ Vector3<real_t> dragForceTang( const Vector3<real_t> & fluidVel, const Vector3<r
 
 // drag correlation based on findings from Felice (1994)
 // used e.g. in Kafui et al (2002)
+inline
 Vector3<real_t> dragForceFelice( const Vector3<real_t> & fluidVel, const Vector3<real_t> & particleVel,
                                  real_t solidVolumeFraction, real_t diameter, real_t fluidDynamicViscosity, real_t fluidDensity )
 {
@@ -180,6 +184,7 @@ Vector3<real_t> dragForceFelice( const Vector3<real_t> & fluidVel, const Vector3
 // drag correlation based on findings from Tenneti, Garg, Subramaniam (2011)
 // used e.g. in Finn, Li, Apte - Particle based modelling and simulation of natural sand dynamics in the wave bottom boundary layer (2016)
 // could be generalized also for non-spherical particles, see Finn et al (2016)
+inline
 Vector3<real_t> dragForceTenneti( const Vector3<real_t> & fluidVel, const Vector3<real_t> & particleVel,
                                   real_t solidVolumeFraction, real_t diameter, real_t fluidDynamicViscosity, real_t fluidDensity )
 {
@@ -211,6 +216,7 @@ Vector3<real_t> dragForceTenneti( const Vector3<real_t> & fluidVel, const Vector
 }
 
 
+inline
 Vector3<real_t> noDragForce( const Vector3<real_t> & /*fluidVel*/, const Vector3<real_t> & /*particleVel*/,
                              real_t /*solidVolumeFraction*/, real_t /*diameter*/, real_t /*fluidDynamicViscosity*/, real_t /*fluidDensity*/ )
 {
diff --git a/src/pe_coupling/discrete_particle_methods/correlations/LiftForceCorrelations.h b/src/pe_coupling/discrete_particle_methods/correlations/LiftForceCorrelations.h
index 8e66b93dab441da4fc7e95bd03f915dac930392a..a5e762d21bfa8ef1b2a08fc68e684ce64f9feeb1 100644
--- a/src/pe_coupling/discrete_particle_methods/correlations/LiftForceCorrelations.h
+++ b/src/pe_coupling/discrete_particle_methods/correlations/LiftForceCorrelations.h
@@ -37,6 +37,7 @@ namespace discrete_particle_methods {
  */
 
 // Saffman lift force
+inline
 Vector3<real_t> liftForceSaffman ( const Vector3<real_t> & fluidVel, const Vector3<real_t> & curlFluidVel, const Vector3<real_t> & particleVel,
                                    real_t diameter, real_t fluidDynamicViscosity, real_t fluidDensity )
 {
@@ -53,7 +54,7 @@ Vector3<real_t> liftForceSaffman ( const Vector3<real_t> & fluidVel, const Vecto
 
 }
 
-Vector3<real_t> noLiftForce ( const Vector3<real_t> &, const Vector3<real_t> &, const Vector3<real_t> &, real_t, real_t, real_t )
+inline Vector3<real_t> noLiftForce ( const Vector3<real_t> &, const Vector3<real_t> &, const Vector3<real_t> &, real_t, real_t, real_t )
 {
    return Vector3<real_t>(real_t(0));
 }
diff --git a/src/pe_coupling/discrete_particle_methods/evaluators/EffectiveViscosityFieldEvaluator.h b/src/pe_coupling/discrete_particle_methods/evaluators/EffectiveViscosityFieldEvaluator.h
index f8039c7e6382121f7e19417476a188ca8cbe8914..10c85b125886e65fabef67be3854560e18fa6060 100644
--- a/src/pe_coupling/discrete_particle_methods/evaluators/EffectiveViscosityFieldEvaluator.h
+++ b/src/pe_coupling/discrete_particle_methods/evaluators/EffectiveViscosityFieldEvaluator.h
@@ -32,20 +32,20 @@ namespace discrete_particle_methods {
 // correlations for the viscosity
 
 // no change in viscosity
-real_t calculateUnchangedEffectiveViscosity( real_t fluidViscosity, real_t /*porosity*/ )
+inline real_t calculateUnchangedEffectiveViscosity( real_t fluidViscosity, real_t /*porosity*/ )
 {
    return fluidViscosity;
 }
 
 // see: Fattahi, Waluga, Wohlmuth - "Large scale lattice Boltzmann simulation for the coupling of free and porous media flow"
-real_t calculateRescaledEffectiveViscosity( real_t fluidViscosity, real_t porosity )
+inline real_t calculateRescaledEffectiveViscosity( real_t fluidViscosity, real_t porosity )
 {
    return fluidViscosity / porosity;
 }
 
 // see: J. R. Finn, M. Li, S. V. Apte - "Particle based modelling and simulation of natural sand dynamics in the wave
 // bottom boundary layer", Journal of Fluid Mechanics 796 (2016) 340–385. doi:10.1017/jfm.2016.246.
-real_t calculateEilersEffectiveViscosity( real_t fluidViscosity, real_t porosity )
+inline real_t calculateEilersEffectiveViscosity( real_t fluidViscosity, real_t porosity )
 {
    const real_t closePackingFraction = real_t(0.64);
    const real_t intrinsicViscosity = real_t(2.5); //for monosized spheres
diff --git a/src/pe_coupling/discrete_particle_methods/evaluators/LubricationForceEvaluator.h b/src/pe_coupling/discrete_particle_methods/evaluators/LubricationForceEvaluator.h
index 2de86148379f3ed7391a09477be7712b44f264e7..7b40b3a57e60e09c5ce7f250874e2d4e5b4f3105 100644
--- a/src/pe_coupling/discrete_particle_methods/evaluators/LubricationForceEvaluator.h
+++ b/src/pe_coupling/discrete_particle_methods/evaluators/LubricationForceEvaluator.h
@@ -88,7 +88,7 @@ private:
 }; // class LubricationForceEvaluator
 
 
-void LubricationForceEvaluator::operator ()()
+inline void LubricationForceEvaluator::operator ()()
 {
    WALBERLA_LOG_PROGRESS( "Calculating Lubrication Force" );
 
@@ -142,7 +142,7 @@ void LubricationForceEvaluator::operator ()()
    }
 }
 
-void LubricationForceEvaluator::treatLubricationSphrSphr( const pe::SphereID sphereI, const pe::SphereID sphereJ, const math::AABB & blockAABB )
+inline void LubricationForceEvaluator::treatLubricationSphrSphr( const pe::SphereID sphereI, const pe::SphereID sphereJ, const math::AABB & blockAABB )
 {
 
    WALBERLA_ASSERT_UNEQUAL( sphereI->getSystemID(), sphereJ->getSystemID() );
@@ -182,7 +182,7 @@ void LubricationForceEvaluator::treatLubricationSphrSphr( const pe::SphereID sph
 
 }
 
-void LubricationForceEvaluator::treatLubricationSphrPlane( const pe::SphereID sphereI, const pe::ConstPlaneID planeJ )
+inline void LubricationForceEvaluator::treatLubricationSphrPlane( const pe::SphereID sphereI, const pe::ConstPlaneID planeJ )
 {
 
    real_t gap = pe::getSurfaceDistance( sphereI, planeJ );
@@ -207,7 +207,7 @@ void LubricationForceEvaluator::treatLubricationSphrPlane( const pe::SphereID sp
 }
 
 
-pe::Vec3 LubricationForceEvaluator::compLubricationSphrSphr( real_t gap, const pe::SphereID sphereI, const pe::SphereID sphereJ) const
+inline pe::Vec3 LubricationForceEvaluator::compLubricationSphrSphr( real_t gap, const pe::SphereID sphereI, const pe::SphereID sphereJ) const
 {
    const pe::Vec3 &posSphereI = sphereI->getPosition();
    const pe::Vec3 &posSphereJ = sphereJ->getPosition();
@@ -260,7 +260,7 @@ pe::Vec3 LubricationForceEvaluator::compLubricationSphrSphr( real_t gap, const p
    return fLub;
 }
 
-pe::Vec3 LubricationForceEvaluator::compLubricationSphrPlane( real_t gap, const pe::SphereID sphereI, const pe::ConstPlaneID planeJ) const
+inline pe::Vec3 LubricationForceEvaluator::compLubricationSphrPlane( real_t gap, const pe::SphereID sphereI, const pe::ConstPlaneID planeJ) const
 {
    const pe::Vec3 &posSphereI( sphereI->getPosition() );
    real_t radiusSphereI = sphereI->getRadius();
diff --git a/src/pe_coupling/geometry/SphereEquivalentDiameter.h b/src/pe_coupling/geometry/SphereEquivalentDiameter.h
index 7e40a843dac7f0d5aca20cb996a5c28eb2914e08..e32dc239388de89c53bbb9b4cffbf7f56352d934 100644
--- a/src/pe_coupling/geometry/SphereEquivalentDiameter.h
+++ b/src/pe_coupling/geometry/SphereEquivalentDiameter.h
@@ -31,7 +31,7 @@ namespace pe_coupling {
 
 
 // calculates sphere-equivalent diameter (diameter of a sphere with same volume as given body)
-real_t getSphereEquivalentDiameter( pe::RigidBody & body )
+inline real_t getSphereEquivalentDiameter( pe::RigidBody & body )
 {
    if( body.getTypeID() == pe::Sphere::getStaticTypeID() || body.getTypeID() == pe::Squirmer::getStaticTypeID() )
    {