Great, please keep me (us) posted. I have several half implemented routines to handle geometric objects (e.g point) with "dimension" and also a set of adapted_boost_array_units.hpp. I am stuck now because I reached the point at which I can't understand the internals of the library.

What I have so far is for example
    using namespace boost::geometry;
    using namespace boost::units;
    model::point<double, 3, cs::units_cartesian<si::length> > p1(1.*si::meter,2.*si::meter,3.*si::meter);
    std::clog << distance(p1,p1) << std::endl; // outputs   0. * meter

and

#include<boost/geometry/geometries/adapted/boost_array_units_cartesian.hpp>
...
    using namespace boost::units;
    boost::array<quantity<si::length>, 3> p5={{10.*si::meter,11.*si::meter,12.*si::meter}};
    std::clog << distance(p5,p5) << std::endl; //outputs 0. * meter

the implementation I have is very primitive and it is practically a mockup because it bypasses most of the features of the Boost.Geometry library; mainly because I didn't understand them internally.

Thank you,
Alfredo