Fix numerous conversions between double and int Vector3D

Code/redblood/Borders.h

@@ -82,26 +82,27 @@ namespace hemelb
         bool wannaSee(value_type const &current) const;
     };
 
-    template<class T> util::Vector3D<T> direction(Borders border)
+    template<class T>
+    util::Vector3D<T> direction(Borders border)
     {
       switch (border)
       {
         case Borders::TOP:
-          return util::Vector3D<T>(1, 0, 0);
+          return {1, 0, 0};
         case Borders::BOTTOM:
-          return util::Vector3D<T>(-1, 0, 0);
+          return {-1, 0, 0};
         case Borders::NORTH:
-          return util::Vector3D<T>(0, 1, 0);
+          return {0, 1, 0};
         case Borders::SOUTH:
-          return util::Vector3D<T>(0, -1, 0);
+          return {0, -1, 0};
         case Borders::WEST:
-          return util::Vector3D<T>(0, 0, -1);
+          return {0, 0, -1};
         case Borders::EAST:
-          return util::Vector3D<T>(0, 0, 1);
+          return {0, 0, 1};
         default:
-          return util::Vector3D<T>(0, 0, 0);
+          return {0, 0, 0};
       };
-      return LatticeVector(0, 0, 0);
+      return {0, 0, 0};
     }
     template<class T> util::Vector3D<T> direction(size_t border)
     {

Code/redblood/Cell.cc

@@ -169,13 +169,13 @@ namespace hemelb
       auto moduli = cell->moduli;
 
       {
-        std::vector<LatticeForceVector> forces(cell->GetNumberOfNodes(), 0.0);
+        std::vector<LatticeForceVector> forces(cell->GetNumberOfNodes(), LatticeForceVector::Zero());
         cell->facetBending(forces);
         push_back_point_data("bending", forces);
       }
 
       {
-        std::vector<LatticeForceVector> forces(cell->GetNumberOfNodes(), 0.0);
+        std::vector<LatticeForceVector> forces(cell->GetNumberOfNodes(), LatticeForceVector::Zero());
         volumeEnergy(cell->GetVertices(),
                      *cell->GetTemplateMesh().GetData(),
                      moduli.volume,
@@ -185,7 +185,7 @@ namespace hemelb
       }
 
       {
-        std::vector<LatticeForceVector> forces(cell->GetNumberOfNodes(), 0.0);
+        std::vector<LatticeForceVector> forces(cell->GetNumberOfNodes(), LatticeForceVector::Zero());
         surfaceEnergy(cell->GetVertices(),
                       *cell->GetTemplateMesh().GetData(),
                       moduli.surface,
@@ -195,7 +195,7 @@ namespace hemelb
       }
 
       {
-        std::vector<LatticeForceVector> forces(cell->GetNumberOfNodes(), 0.0);
+        std::vector<LatticeForceVector> forces(cell->GetNumberOfNodes(), LatticeForceVector::Zero());
         strainEnergy(cell->GetVertices(),
                      *cell->GetTemplateMesh().GetData(),
                      moduli.strain,

Code/redblood/CellArmy.h

@@ -338,7 +338,8 @@ namespace hemelb
       auto const barycenter = cell->GetBarycenter();
 
       //! @todo: #623 AddCell should only be called if the subdomain contains the relevant RBC inlet
-      auto const iter = globalCoordsToProcMap.find(barycenter);
+      // TODO: #759 truncation of barycenter
+      auto const iter = globalCoordsToProcMap.find(LatticeVector{barycenter});
       bool insertAtThisRank = (iter != globalCoordsToProcMap.end()) && (iter->second == neighbourDependenciesGraph.Rank());
       if (insertAtThisRank)
       {

Code/redblood/GridAndCell.h

@@ -69,7 +69,7 @@ namespace hemelb
     Dimensionless forcesOnGrid(std::shared_ptr<CellBase const> cell,
                                geometry::LatticeData &latticeData)
     {
-      std::vector<LatticeForceVector> forces(cell->GetNumberOfNodes(), 0);
+      std::vector<LatticeForceVector> forces(cell->GetNumberOfNodes(), {0.0, 0.0, 0.0});
       return forcesOnGrid<LATTICE, STENCIL>(cell, forces, latticeData);
     }
   }

Code/redblood/VertexBag.cc

@@ -13,7 +13,7 @@ namespace hemelb
     {
       // avoids a warning
 #     ifndef HEMELB_DOING_UNITTESTS
-      proc_t get_proc(geometry::LatticeData const &latDat, LatticePosition const &pos)
+      proc_t get_proc(geometry::LatticeData const &latDat, LatticeVector const &pos)
       {
         proc_t procid;
         site_t siteid;
@@ -104,7 +104,7 @@ namespace hemelb
         std::shared_ptr<CellBase const> cell, geometry::LatticeData const &latticeData,
         proc_t selfRegion)
     {
-      auto proc_getter = [&latticeData](LatticePosition const &position)
+      auto proc_getter = [&latticeData](LatticeVector const &position)
       {
         return get_proc(latticeData, position);
       };

Code/unittests/FourCubeLatticeData.h

@@ -108,24 +108,25 @@ namespace hemelb
 
                   float randomDistance = (float(std::rand() % 10000) / 10000.0);
 
+		  using CutType = io::formats::geometry::CutType;
                   // The inlet is by the minimal z value.
                   if (neighK < minInd)
                   {
                     link.ioletId = 0;
-                    link.type = geometry::GeometrySiteLink::INLET_INTERSECTION;
+                    link.type = CutType::INLET;
                     link.distanceToIntersection = randomDistance;
                   }
                   // The outlet is by the maximal z value.
                   else if (neighK > maxInd)
                   {
                     link.ioletId = 0;
-                    link.type = geometry::GeometrySiteLink::OUTLET_INTERSECTION;
+                    link.type = CutType::OUTLET;
                     link.distanceToIntersection = randomDistance;
                   }
                   // Walls are by extremes of x and y.
                   else if (neighI < minInd || neighJ < minInd || neighI > maxInd || neighJ > maxInd)
                   {
-                    link.type = geometry::GeometrySiteLink::WALL_INTERSECTION;
+                    link.type = CutType::WALL;
                     link.distanceToIntersection = randomDistance;
                   }
 

Code/unittests/helpers/LatticeDataAccess.h

@@ -90,7 +90,7 @@ namespace hemelb
       }
       inline void LatticeDataAccess::ZeroOutForces() const
       {
-        std::fill(latDat->forceAtSite.begin(), latDat->forceAtSite.end(), 0);
+        std::fill(latDat->forceAtSite.begin(), latDat->forceAtSite.end(), LatticeForceVector::Zero());
       }
 
       template<class LATTICE>

Code/unittests/redblood/BendingTests.h

@@ -55,7 +55,7 @@ namespace hemelb
           {
             for (auto const theta : { 1e-2, 2e-2, 3e-2 })
             {
-              std::fill(forces.begin(), forces.end(), 0e0);
+              std::fill(forces.begin(), forces.end(), LatticeForceVector::Zero());
               vertices.back() = bending(mesh.vertices.back(), theta);
               CPPUNIT_ASSERT_DOUBLES_EQUAL(std::sqrt(3.) * theta * theta * moduli, energy(), 1e-10);
               CPPUNIT_ASSERT_DOUBLES_EQUAL(std::sqrt(3.) * theta * theta * moduli, energy(), 1e-10);
@@ -166,7 +166,7 @@ namespace hemelb
             auto const theta = orientedAngle(Facet(vertices, mesh.facets, 0),
                                              Facet(vertices, mesh.facets, 1));
             auto const epsilon = 1e-4;
-            std::fill(forces.begin(), forces.end(), 0e0);
+            std::fill(forces.begin(), forces.end(), LatticeForceVector::Zero());
             auto const oldPos = vertices[node];
 
             auto const e0 = energyAndForces();

Code/unittests/redblood/CellCellInteractionTests.h

@@ -141,13 +141,13 @@ namespace hemelb
         CPPUNIT_ASSERT(std::distance(empty.first, empty.second) == 0);
         CPPUNIT_ASSERT(std::distance(haloed.first, haloed.second) == 1);
 
-        CPPUNIT_ASSERT(helpers::is_zero(omega.first->first));
+        CPPUNIT_ASSERT_EQUAL(LatticeVector::Zero(), omega.first->first);
         CPPUNIT_ASSERT(omega.first->second.cellIterator == cells.begin());
         CPPUNIT_ASSERT(omega.first->second.nodeIndex == 0 or omega.first->second.nodeIndex == 1);
         CPPUNIT_ASSERT(omega.first->second.nearBorder == (size_t) Borders::CENTER);
         DivideConquer<CellReference>::const_iterator other = omega.first;
         ++other;
-        CPPUNIT_ASSERT(helpers::is_zero(other->first));
+        CPPUNIT_ASSERT_EQUAL(LatticeVector::Zero(), other->first);
         CPPUNIT_ASSERT(other->second.cellIterator == cells.begin());
         CPPUNIT_ASSERT(other->second.nodeIndex == 0 or other->second.nodeIndex == 1);
         CPPUNIT_ASSERT(other->second.nodeIndex != omega.first->second.nodeIndex);

Code/unittests/redblood/CellForceSpreadTests.h

@@ -76,7 +76,8 @@ namespace hemelb
                                                                                             {}),
                                                                   details::SpreadForces(forces,
                                                                                         *latDat));
-        return latDat->GetSite(center).GetForce();
+	// TODO #759: is this truncation OK?
+        return latDat->GetSite(LatticeVector{center}).GetForce();
       }
 
       void CellForceSpreadTests::setUp()

Code/unittests/redblood/CellInserterTests.h

@@ -35,7 +35,7 @@ namespace hemelb
 
           void setUp()
           {
-            converter.reset(new util::UnitConverter(0.5, 0.6, 0e0, 1000.0, 0.0));
+            converter.reset(new util::UnitConverter(0.5, 0.6, PhysicalPosition::Zero(), 1000.0, 0.0));
             every = 10;
             offset = 5;
             cells.emplace("joe", std::make_shared<Cell>(tetrahedron()));

Code/unittests/redblood/CellTests.h

@@ -51,6 +51,11 @@ namespace hemelb
               }
           };
 
+	template <typename VEC3>
+	auto vector_of_zero_vec3(std::size_t count) {
+	  return std::vector<VEC3>{count, VEC3::Zero()};
+	}
+
         public:
           void setUp()
           {
@@ -69,7 +74,7 @@ namespace hemelb
           void testNullTemplateScaling()
           {
             Mesh templateMesh(original);
-            std::vector<LatticeForceVector> forces(original.vertices.size(), 0);
+	    auto forces = vector_of_zero_vec3<LatticeForceVector>(original.vertices.size());
 
             std::vector<Dimensionless> scales;
             scales.push_back(1.0);
@@ -100,8 +105,8 @@ namespace hemelb
             scaled.GetVertices()[1] += LatticePosition(0.0837, -0.012632, 0.0872935);
             scaled.GetVertices()[2] += LatticePosition(0.02631, -0.00824223, -0.098362);
 
-            std::vector<LatticeForceVector> uforces(original.vertices.size(), 0),
-                sforces(original.vertices.size(), 0);
+            auto uforces = vector_of_zero_vec3<LatticeForceVector>(original.vertices.size());
+	    auto sforces = vector_of_zero_vec3<LatticeForceVector>(original.vertices.size());
 
             scaled *= 1.1;
             auto const uenergy = GetCellWithEnergy(scaled, templateMesh, 1.1)(uforces);

Code/unittests/redblood/DivideConquerTests.h

@@ -59,7 +59,7 @@ namespace hemelb
                                             LatticeVector(1, -2, 2) };
 
         for (size_t i(0); i < N; ++i)
-          CPPUNIT_ASSERT(helpers::is_zero(dnc.DowngradeKey(inputs[i]) - expected[i]));
+          CPPUNIT_ASSERT_EQUAL(LatticeVector::Zero(), dnc.DowngradeKey(inputs[i]) - expected[i]);
       }
 
       void DivideAndConquerTests::testNoDowngradeKey()

Code/unittests/redblood/Fixtures.h

@@ -105,7 +105,7 @@ namespace hemelb
           FourCubeBasedTestFixture::setUp();
           mesh = redblood::refine(initial_mesh(), refinement());
           mesh *= Dimensionless(CubeSize() - 3) * 0.5;
-          mesh += Dimensionless(CubeSize()) * 0.5;
+          mesh += LatticePosition{Dimensionless(CubeSize()) * 0.5};
         }
 
       protected:

Code/unittests/redblood/FlowExtensionTests.h

@@ -132,7 +132,7 @@ namespace hemelb
                       "    <fadelength units=\"m\" value=\"0.05\" />"
                       "  </flowextension>"
                       "</inlet>");
-            util::UnitConverter converter(0.5, 0.6, 0.7, DEFAULT_FLUID_DENSITY_Kg_per_m3, 0.0);
+            util::UnitConverter converter(0.5, 0.6, PhysicalPosition{0.7}, DEFAULT_FLUID_DENSITY_Kg_per_m3, 0.0);
             auto flow = readFlowExtension(doc.FirstChildElement("inlet"), converter);
             // Normal is opposite direction compared to XML inlet definition
             CPPUNIT_ASSERT_DOUBLES_EQUAL(-2e0 / std::sqrt(2),

Code/unittests/redblood/InterpolationTests.h

@@ -56,6 +56,11 @@ namespace hemelb
                                        LatticePosition(-1, 2, 1).Dot(pos),
                                        LatticePosition(0, 0, 1).Dot(pos));
               }
+              // We need to accept grid point and continuous input points
+              LatticePosition operator()(LatticeVector const &pos) const
+              {
+                return operator()(pos.as<double>());
+              }
               LatticePosition operator()(Dimensionless x, Dimensionless y, Dimensionless z) const
               {
                 return operator()(LatticePosition(x, y, z));
@@ -64,6 +69,7 @@ namespace hemelb
 
           struct QuadraticFunction
           {
+              // We need to accept grid point and continuous input points
               LatticePosition operator()(LatticePosition const &pos) const
               {
                 Dimensionless const offset(0);
@@ -74,6 +80,11 @@ namespace hemelb
                                        (LatticePosition(0, 0, 1).Dot(pos) - offset)
                                            * (LatticePosition(0, 0, 1).Dot(pos) - offset));
               }
+              // We need to accept grid point and continuous input points
+              LatticePosition operator()(LatticeVector const &pos) const
+              {
+                return operator()(pos.as<double>());
+              }
               LatticePosition operator()(Dimensionless x, Dimensionless y, Dimensionless z) const
               {
                 return operator()(LatticePosition(x, y, z));
@@ -102,7 +113,7 @@ namespace hemelb
               for (size_t j(0); j < incs[i]; ++j, ++iterator)
                 ;
 
-              CPPUNIT_ASSERT(helpers::is_zero(*iterator - vectors[i]));
+              CPPUNIT_ASSERT_EQUAL(LatticeVector::Zero(), (*iterator - vectors[i]));
               CPPUNIT_ASSERT(iterator.IsValid());
             }
 

Code/unittests/redblood/VertexBagTests.h

@@ -114,7 +114,7 @@ namespace hemelb
                                  splits[0]->GetVertices().size() + splits[2]->GetVertices().size());
           }
 
-          static proc_t proc_at_pos(LatticePosition const &position)
+          static proc_t proc_at_pos(LatticeVector const &position)
           {
             return (position.x > 0 ?
               0 :

Code/unittests/redblood/parallel/Fixtures.h

@@ -19,7 +19,7 @@ namespace hemelb
       //! \brief gathers mid-domain and egde positions from all procs
       //! \details If there are insufficient number of edges, mid-domains are used instead.
       //! erase removes the components from the first process.
-      std::vector<LatticePosition> GatherSpecialPositions(geometry::LatticeData const & latDat,
+      std::vector<LatticeVector> GatherSpecialPositions(geometry::LatticeData const & latDat,
                                                           size_t mid, size_t edges,
                                                           net::MpiCommunicator const &c)
       {
@@ -28,7 +28,7 @@ namespace hemelb
 
         int const nMids = latDat.GetMidDomainCollisionCount(0);
         int const nEdges = latDat.GetDomainEdgeCollisionCount(0);
-        std::vector<LatticePosition> positions(c.Size() * (mid + edges));
+        std::vector<LatticeVector> positions(c.Size() * (mid + edges));
         std::vector<int> shuf(nMids);
         std::iota(shuf.begin(), shuf.end(), 0);
         std::shuffle(shuf.begin(), shuf.end(), g);
@@ -106,7 +106,7 @@ namespace hemelb
           LatticeEnergy operator()(std::vector<LatticeForceVector> &f) const override
           {
             f.resize(GetNumberOfNodes());
-            std::fill(f.begin(), f.end(), force);
+            std::fill(f.begin(), f.end(), LatticeForceVector{force});
             return 0e0;
           }
           std::unique_ptr<CellBase> cloneImpl() const override

Code/unittests/redblood/parallel/MPIIntegrateVelocities.h

@@ -187,7 +187,7 @@ namespace hemelb
           typedef std::vector<LatticePosition> Positions;
           Positions positions = world.Rank() == 0 ?
             cell->GetVertices() :
-            Positions(cell->GetNumberOfNodes(), 0);
+            Positions(cell->GetNumberOfNodes(), LatticePosition::Zero());
           world.Broadcast(positions, 0);
           if (world.Rank() != 0 and owned.find(cell) != owned.end())
           {

Code/unittests/redblood/parallel/MPILockStepTests.h

@@ -140,7 +140,7 @@ namespace hemelb
           geometry::LatticeData const &latDat,
           std::size_t &nbtests,
           hemelb::redblood::CellContainer const &cells) {
-        std::vector<LatticePosition> positions;
+        std::vector<LatticeVector> positions;
         std::vector<LatticeForceVector> forces;
         for(site_t i(0); i < latDat.GetLocalFluidSiteCount(); ++i)
         {
@@ -157,7 +157,7 @@ namespace hemelb
         }
         if(world.Rank() == 0)
         {
-          auto const parallel_positions = world.Gather(std::vector<LatticePosition>{}, 0);
+          auto const parallel_positions = world.Gather(std::vector<LatticeVector>{}, 0);
           auto const parallel_forces = world.Gather(std::vector<LatticeForceVector>{}, 0);
           CPPUNIT_ASSERT_EQUAL(positions.size(), parallel_positions.size());
           for(std::size_t i(0); i < positions.size(); ++i)

Code/unittests/redblood/parallel/MPISpreadForcesTests.h

@@ -176,7 +176,7 @@ namespace hemelb
             }
             else
             {
-              forces.push_back(0);
+              forces.push_back(LatticeForceVector::Zero());
             }
           }
         }

Code/unittests/redblood/parallel/ParallelFixtureTests.h

@@ -118,7 +118,7 @@ namespace hemelb
         for(std::size_t i(0); i < positions.size(); ++i)
         {
           auto const procs = hemelb::redblood::parallel::details::positionAffectsProcs<Stencil>(
-              globalCoordsToProcMap, positions[i]);
+												globalCoordsToProcMap, positions[i].as<LatticeDistance>());
 
           // Send set of affected procs as known by owner proc
           decltype(world.Rank()) positions_are_from_this_proc = i / (nmid + nedges) + 1;