diff --git a/include/diy/decomposition.hpp b/include/diy/decomposition.hpp
index f719f620452cb6341d9b68beff58701c6045ee04..7f97e867442c1a3d67694a9923c8a6b3a5291cbf 100644
--- a/include/diy/decomposition.hpp
+++ b/include/diy/decomposition.hpp
@@ -123,6 +123,7 @@ namespace detail
     template<class Point>
     int             lowest_gid(const Point& p) const;
 
+    DivisionsVector gid_to_coords(int gid) const                                { DivisionsVector coords; gid_to_coords(gid, coords); return coords; }
     void            gid_to_coords(int gid, DivisionsVector& coords) const       { gid_to_coords(gid, coords, divisions); }
     int             coords_to_gid(const DivisionsVector& coords) const          { return coords_to_gid(coords, divisions); }
     void            fill_divisions(std::vector<int>& divisions) const;
@@ -411,6 +412,7 @@ void
 diy::RegularDecomposer<Bounds>::
 gid_to_coords(int gid, DivisionsVector& coords, const DivisionsVector& divisions)
 {
+  coords.clear();
   int dim = static_cast<int>(divisions.size());
   for (int i = 0; i < dim; ++i)
   {
diff --git a/tests/decomposer.cpp b/tests/decomposer.cpp
index e8f928199d368da906037c636dd60b3d09d24bc6..a0529a06396b5c6ecb7201ff6724f6209d4772ac 100644
--- a/tests/decomposer.cpp
+++ b/tests/decomposer.cpp
@@ -46,6 +46,14 @@ TEST_CASE("RegularDecomposer", "[decomposition]")
 
         diy::ContiguousAssigner assigner(nblocks, nblocks);
         decomposer.decompose(0, assigner, test);
+
+        // test the return-by-value gid_to_coords
+        auto coords = decomposer.gid_to_coords(2);
+        REQUIRE(coords == diy::RegularDecomposer<diy::ContinuousBounds>::DivisionsVector { 2, 0, 0 });
+
+        // test that gid_to_coords clears the vector
+        decomposer.gid_to_coords(3, coords);
+        REQUIRE(coords == diy::RegularDecomposer<diy::ContinuousBounds>::DivisionsVector { 0, 1, 0 });
     }
 
     SECTION("1D decomposition of an interval")