more functionality for geometry

Framework/Geometry/QuantityVector.h

@@ -11,6 +11,7 @@ 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));
     
 public:
     Eigen::Vector3d eVector;
@@ -25,16 +26,21 @@ public:
     {
     }
     
-    Quantity operator[](size_t index) const
+    auto operator[](size_t index) const
     {
         return Quantity(phys::units::detail::magnitude_tag, eVector[index]);
     }
     
-    Quantity norm() const
+    auto norm() const
     {
         return Quantity(phys::units::detail::magnitude_tag, eVector.norm());
     }
     
+    auto squaredNorm() const
+    {
+        return QuantitySquared(phys::units::detail::magnitude_tag, eVector.squaredNorm());
+    }
+    
     auto operator+(QuantityVector<dim> const& pQVec) const
     {
         return QuantityVector<dim>(eVector + pQVec.eVector);

Framework/Geometry/Sphere.h

+#ifndef SPHERE_H_
+#define SPHERE_H_
+
+#include "Point.h"
+#include <phys/units/quantity.hpp>
+
+class Sphere
+{
+    using Length = phys::units::quantity<phys::units::length_d, double>;
+    Point center;
+    Length const radius;
+    
+public:
+    Sphere(Point const& pCenter, Length const pRadius) :
+        center(pCenter), radius(pRadius)
+    {
+        
+    }
+        
+    auto isInside(Point const& p) const
+    {
+        return radius*radius > (center - p).squaredNorm();
+    }
+
+};
+
+#endif

Framework/Geometry/Vector.h

@@ -49,6 +49,28 @@ public:
         return BaseVector<dim>::qVector.norm();
     }
     
+    auto squaredNorm() const
+    {
+        return BaseVector<dim>::qVector.squaredNorm();
+    }
+        
+    template <typename dim2>
+    auto parallelProjectionOnto(BaseVector<dim2> const& pVec, CoordinateSystem const& pCS) const
+    {
+        auto const ourCompVec = getComponents(pCS);
+        auto const otherCompVec = pVec.getComponents(pVec);
+        auto const& a = ourCompVec.eVector;
+        auto const& b = otherCompVec.eVector;
+        
+        return Vector<dim>(pCS, (a * b) / b.squaredNorm() * b);
+    }
+    
+    template <typename dim2>
+    auto parallelProjectionOnto(BaseVector<dim2> const& pVec)
+    {
+        return parallelProjectionOnto<dim2>(pVec, *BaseVector<dim>::cs);
+    }
+    
     //~ template <typename dim2>
     //~ auto operator*(Vector<dim2> const& pVec)
     //~ {