diff --git a/Framework/Geometry/QuantityVector.h b/Framework/Geometry/QuantityVector.h
index 64a312f3ddbabc4afee627f87b0f09300055bbd0..36cf0d567f520f8ba93f51cbbf73df9b92ca59a3 100644
--- a/Framework/Geometry/QuantityVector.h
+++ b/Framework/Geometry/QuantityVector.h
@@ -5,16 +5,19 @@
 
 #include <Eigen/Dense>
 #include <iostream>
+#include <utility>
 
 template <typename dim>
 class QuantityVector
 {
 protected:
     using Quantity = phys::units::quantity<dim, double>;
-    using QuantitySquared = decltype(Quantity(phys::units::detail::magnitude_tag, 0) * Quantity(phys::units::detail::magnitude_tag, 0));
+    //using QuantitySquared = decltype(std::declval<Quantity>() * std::declval<Quantity>());
     
 public:
     Eigen::Vector3d eVector;
+    
+    typedef dim dimension;
 
     QuantityVector(Quantity a, Quantity b, Quantity c) :
         eVector{a.magnitude(), b.magnitude(), c.magnitude()}
@@ -38,6 +41,7 @@ public:
     
     auto squaredNorm() const
     {
+        using QuantitySquared = decltype(std::declval<Quantity>() * std::declval<Quantity>());
         return QuantitySquared(phys::units::detail::magnitude_tag, eVector.squaredNorm());
     }
     
@@ -51,14 +55,56 @@ public:
         return QuantityVector<dim>(eVector - pQVec.eVector);
     }
     
-    //auto operator*(
+    template <typename ScalarDim>
+    auto operator*(phys::units::quantity<ScalarDim, double> const p) const
+    {
+        return QuantityVector<typename phys::units::detail::Product<ScalarDim, dim, double, double>::dimension_type>(eVector * p.magnitude());
+        // TODO: this function does not work if the result is dimensionless, as
+        // dimensionless quantities are "cast" back to plain old double in PhysUnits.
+        // Either change PhysUnits, or cover this case with a template specialization?
+    }
+    
+    auto operator*(double const p) const
+    {
+        return QuantityVector<dim>(eVector * p);
+    }
+    
+    auto& operator*=(double const p)
+    {
+        eVector *= p;
+        return *this;
+    }
+    
+    auto& operator+=(QuantityVector<dim> const& pQVec)
+    {
+        eVector += pQVec.eVector;
+        return *this;
+    }
+    
+    auto& operator-=(QuantityVector<dim> const& pQVec)
+    {
+        eVector -= pQVec.eVector;
+        return *this;
+    }
+    
+    auto& operator-() const
+    {
+        return QuantityVector<dim>(-eVector);
+    }
+    
+    auto normalized() const
+    {
+        return (*this) * (1 / norm());
+    }
 };
 
 template <typename dim>
 auto& operator<<(std::ostream& os, QuantityVector<dim> qv)
 {
+    using Quantity = phys::units::quantity<dim, double>;
+    
     os << '(' << qv.eVector(0) << ' ' << qv.eVector(1) << ' ' << qv.eVector(2)
-       << ") " << phys::units::to_unit_symbol(phys::units::quantity<dim, double>());
+       << ") " << phys::units::to_unit_symbol<dim, double>(Quantity(phys::units::detail::magnitude_tag, 1));
     return os;
 }
 
diff --git a/Framework/Geometry/Vector.h b/Framework/Geometry/Vector.h
index 417cb19a52cc5bae2af85c488ad09db766f7dd2d..47b91f7ccc6ce6b3075a151f9b7085abf70e1b59 100644
--- a/Framework/Geometry/Vector.h
+++ b/Framework/Geometry/Vector.h
@@ -71,6 +71,58 @@ public:
         return parallelProjectionOnto<dim2>(pVec, *BaseVector<dim>::cs);
     }
     
+    auto operator+(Vector<dim> const& pVec) const
+    {
+        auto const components = getComponents(*BaseVector<dim>::cs) + pVec.getComponents(*BaseVector<dim>::cs);
+        return Vector<dim>(*BaseVector<dim>::cs, components);
+    }
+    
+    auto operator-(Vector<dim> const& pVec) const
+    {
+        auto const components = getComponents() - pVec.getComponents(*BaseVector<dim>::cs);
+        return Vector<dim>(*BaseVector<dim>::cs, components);
+    }
+    
+    auto& operator*=(double const p)
+    {
+        BaseVector<dim>::qVector *= p;
+        return *this;
+    }
+    
+    template <typename ScalarDim>
+    auto operator*(phys::units::quantity<ScalarDim, double> const p) const
+    {
+        using res_dim = typename decltype(BaseVector<dim>::qVector * p)::dimension;
+        return Vector<res_dim>(*BaseVector<dim>::cs, BaseVector<dim>::qVector * p);        
+    }
+    
+    auto operator*(double const p) const
+    {        
+        return Vector<dim>(*BaseVector<dim>::cs, BaseVector<dim>::qVector * p);
+    }
+    
+    auto& operator+=(Vector<dim> const& pVec)
+    {
+        BaseVector<dim>::qVector += pVec.getComponents(*BaseVector<dim>::cs);
+        return *this;
+    }
+    
+    auto& operator-=(Vector<dim> const& pVec)
+    {
+        BaseVector<dim>::qVector -= pVec.getComponents(*BaseVector<dim>::cs);
+        return *this;
+    }
+    
+    auto& operator-() const
+    {
+        return Vector<dim>(*BaseVector<dim>::cs, - BaseVector<dim>::qVector);
+    }
+    
+    auto normalized() const
+    {
+        return (*this) * (1 / norm());
+    }
+    
     //~ template <typename dim2>
     //~ auto operator*(Vector<dim2> const& pVec)
     //~ {