Somewhat better wrapped support for Icosahedral

First cut at stencil

Grid/cartesian/CartesianCrossIcosahedron.h

@@ -42,7 +42,10 @@ enum NorthSouth {
   South = 0
 };
 
-const int num_icosahedron_tiles = 10;
+const int IcosahedralPatches = 10;
+const int HemiPatches=IcosahedralPatches/2;
+const int NorthernHemisphere = HemiPatches;
+const int SouthernHemisphere = 0;
 
 class GridCartesianCrossIcosahedron: public GridCartesian {
 
@@ -86,9 +89,12 @@ public:
 
     assert(simd_layout[0]==1); // Force simd into perpendicular dimensions
     assert(simd_layout[1]==1); // to avoid pole storage complexity interacting with SIMD.
-    assert(dimensions[_ndimension-1]==num_icosahedron_tiles);
+    assert(dimensions[_ndimension-1]==IcosahedralPatches);
     assert(processor_grid[_ndimension-1]<=2); // Keeps the patches that need a pole on the same node
 
+    // Save a copy of the basic cartesian initialisation volume
+    cartesianOsites = this->_osites;
+
     // allocate the pole storage if we are seeking vertex domain data
     if ( meshType == IcosahedralVertices ) {
       InitPoles();
@@ -106,16 +112,18 @@ public:
   int southPoleOsite;
   int northPoleOsites;
   int southPoleOsites;
+  int cartesianOsites;
 
   virtual int isIcosahedral(void)           override { return 1;}
   virtual int isIcosahedralVertex(void)     override { return meshType==IcosahedralVertices;}
   virtual int isIcosahedralEdge  (void)     override { return meshType==IcosahedralEdges;}
-  virtual int ownsNorthPole(void)   const override { return hasNorthPole; };
   virtual int NorthPoleOsite(void)  const override { return northPoleOsite; };
   virtual int NorthPoleOsites(void) const override { return northPoleOsites; };
-  virtual int ownsSouthPole(void)   const override { return hasSouthPole; };
   virtual int SouthPoleOsite(void)  const override { return southPoleOsite; };
   virtual int SouthPoleOsites(void) const override { return southPoleOsites; };
+  virtual int ownsNorthPole(void)   const override { return hasNorthPole; };
+  virtual int ownsSouthPole(void)   const override { return hasSouthPole; };
+  virtual int CartesianOsites(void) const override { return cartesianOsites; };
 
   void InitPoles(void)
   {
@@ -166,7 +174,7 @@ public:
      * Hence all 5 patches associated with the pole must have the
      * appropriate "corner" of the patch L^2 located on the SAME rank.
      */ 
-     
+    
     if( (pcoor[xdim]==pgrid[xdim]-1) && (pcoor[ydim]==0) && (pcoor[Ndm1]==0) ){
       hasSouthPole   =1;
       southPoleOsite=this->_osites;

Grid/cartesian/Cartesian_base.h

@@ -96,6 +96,7 @@ public:
   virtual int SouthPoleOsite(void) const { return 0; };
   virtual int NorthPoleOsites(void) const { std::cout << "base osites" <<std::endl;return 0; };
   virtual int SouthPoleOsites(void) const { std::cout << "base osites" <<std::endl;return 0; };
+  virtual int CartesianOsites(void) const { return this->oSites(); };
 
   ////////////////////////////////////////////////////////////////
   // Checkerboarding interface is virtual and overridden by 

Grid/lattice/Lattice_base.h

@@ -391,7 +391,7 @@ template<class vobj> std::ostream& operator<< (std::ostream& stream, const Latti
     stream<<"]\t";
     stream<<ss<<std::endl;
   }
-  if ( o.Grid()->isIcosahedral() ) {
+  if ( o.Grid()->isIcosahedralVertex() ) {
     uint64_t psites=1;
     Coordinate perpdims;
     for(int d=2;d<o.Grid()->_ndimension-1;d++){

Grid/lattice/Lattice_peekpoke.h

@@ -281,6 +281,131 @@ void pokePole(const sobj &s,Lattice<vobj> &l,const Coordinate &orthog,NorthSouth
   return;
 };
 
+
+template<class vobj,class sobj>
+void peekLocalPole(sobj &s,const Lattice<vobj> &l,const Coordinate &orthog,NorthSouth isNorth)
+{
+  s=Zero();
+  
+  GridBase *grid=l.Grid();
+
+  assert(grid->isIcosahedral());
+  assert(grid->isIcosahedralVertex());
+
+  int Nsimd = grid->Nsimd();
+
+  int rank;
+
+  int Ndm1         = grid->_ndimension-1;
+  Coordinate pgrid = grid->ProcessorGrid();
+  const int xdim=0;
+  const int ydim=1;
+  const int pdim=Ndm1;
+
+  int64_t pole_osite;
+  int64_t pole_isite;
+  Coordinate rdims;
+  Coordinate idims;
+  Coordinate ocoor;
+  Coordinate icoor;
+  //  Coordinate pcoor(grid->_ndimension);
+  for(int d=2;d<Ndm1;d++){
+    int dd=d-2;
+    rdims.push_back(grid->_rdimensions[d]);
+    idims.push_back(grid->_simd_layout[d]);
+    icoor.push_back((orthog[dd]%grid->_ldimensions[d])/grid->_rdimensions[d]);
+    ocoor.push_back(orthog[dd]%grid->_rdimensions[d]);
+    //    pcoor[d] = orthog[dd]/grid->_ldimensions[d];
+  }
+  Lexicographic::IndexFromCoor(ocoor,pole_osite,rdims);
+  Lexicographic::IndexFromCoor(icoor,pole_isite,idims);
+  
+  int64_t osite;
+  if(isNorth == North){
+    //    pcoor[xdim] = 0;
+    //    pcoor[ydim] = pgrid[ydim]-1;
+    //    pcoor[Ndm1] = pgrid[Ndm1]-1;
+    osite = pole_osite + grid->NorthPoleOsite();
+    assert(grid->ownsNorthPole());
+  } else {
+    //    pcoor[xdim] = pgrid[xdim]-1;
+    //    pcoor[ydim] = 0;
+    //    pcoor[Ndm1] = 0;
+    osite = pole_osite + grid->SouthPoleOsite();
+    assert(grid->ownsSouthPole());
+  }
+
+  ExtractBuffer<sobj> buf(Nsimd);
+  autoView( l_v , l, CpuWrite);
+  extract(l_v[osite],buf);
+  s = buf[pole_isite];
+
+  return;
+};
+template<class vobj,class sobj>
+void pokeLocalPole(const sobj &s,Lattice<vobj> &l,const Coordinate &orthog,NorthSouth isNorth)
+{
+  GridBase *grid=l.Grid();
+
+  assert(grid->isIcosahedral());
+  assert(grid->isIcosahedralVertex());
+
+  int Nsimd = grid->Nsimd();
+  int rank;
+  int Ndm1         = grid->_ndimension-1;
+
+  const int xdim=0;
+  const int ydim=1;
+  const int pdim=Ndm1;
+
+  int64_t pole_osite;
+  int64_t pole_isite;
+  Coordinate rdims;
+  Coordinate idims;
+  Coordinate ocoor;
+  Coordinate icoor;
+  //  Coordinate pcoor(grid->_ndimension,0);
+  for(int d=2;d<Ndm1;d++){
+    int dd = d-2;
+    rdims.push_back(grid->_rdimensions[d]);
+    idims.push_back(grid->_simd_layout[d]);
+    icoor.push_back((orthog[dd]%grid->_ldimensions[d])/grid->_rdimensions[d]);
+    ocoor.push_back(orthog[dd]%grid->_rdimensions[d]);
+    //    pcoor[d] = orthog[dd]/grid->_ldimensions[d];
+
+    int o = orthog[dd];
+    int r = grid->_rdimensions[d];
+    int omr = o % r;
+  }
+  Lexicographic::IndexFromCoor(ocoor,pole_osite,rdims);
+  Lexicographic::IndexFromCoor(icoor,pole_isite,idims);
+  
+  int64_t osite;
+  int insert=0;
+  if(isNorth ==North){
+    //    pcoor[xdim] = 0;
+    //    pcoor[ydim] = pgrid[ydim]-1;
+    //    pcoor[Ndm1] = pgrid[Ndm1]-1;
+    osite = pole_osite + grid->NorthPoleOsite();
+    assert(grid->ownsNorthPole());
+  } else {
+    //    pcoor[xdim] = pgrid[xdim]-1;
+    //    pcoor[ydim] = 0;
+    //    pcoor[Ndm1] = 0;
+    osite = pole_osite + grid->SouthPoleOsite();
+    assert(grid->ownsSouthPole());
+  }
+
+  // extract-modify-merge cycle is easiest way and this is not perf critical
+  ExtractBuffer<sobj> buf(Nsimd);
+  autoView( l_v , l, CpuWrite);
+  extract(l_v[osite],buf);
+  buf[pole_isite] = s;
+  merge(l_v[osite],buf);
+  
+  return;
+};
+
 //////////////////////////////////////////////////////////
 // Peek a scalar object from the SIMD array
 //////////////////////////////////////////////////////////

Grid/lattice/Lattice_rng.h

@@ -48,31 +48,45 @@ NAMESPACE_BEGIN(Grid);
 //////////////////////////////////////////////////////////////
 inline int RNGfillable(GridBase *coarse,GridBase *fine)
 {
+  if ( coarse == fine ) return 1;
 
-  int rngdims = coarse->_ndimension;
-
-  // trivially extended in higher dims, with locality guaranteeing RNG state is local to node
-  int lowerdims   = fine->_ndimension - coarse->_ndimension;
-  assert(lowerdims >= 0);
-  for(int d=0;d<lowerdims;d++){
-    assert(fine->_simd_layout[d]==1);
-    assert(fine->_processors[d]==1);
+  if ( coarse->isIcosahedral()) assert(coarse->isIcosahedralEdge());
+  
+  if ( fine->isIcosahedralVertex() && coarse->isIcosahedralEdge() ) {
+    assert(fine->Nd()==coarse->Nd());
+    for(int d=0;d<fine->Nd();d++){
+      assert(fine->LocalDimensions()[d] == coarse->LocalDimensions()[d]);
+    }
+    return 1;
   }
+    
+  {
+    
+    int rngdims = coarse->_ndimension;
 
-  int multiplicity=1;
-  for(int d=0;d<lowerdims;d++){
-    multiplicity=multiplicity*fine->_rdimensions[d];
-  }
-  // local and global volumes subdivide cleanly after SIMDization
-  for(int d=0;d<rngdims;d++){
-    int fd= d+lowerdims;
-    assert(coarse->_processors[d]  == fine->_processors[fd]);
-    assert(coarse->_simd_layout[d] == fine->_simd_layout[fd]);
-    assert(((fine->_rdimensions[fd] / coarse->_rdimensions[d])* coarse->_rdimensions[d])==fine->_rdimensions[fd]); 
+    // trivially extended in higher dims, with locality guaranteeing RNG state is local to node
+    int lowerdims   = fine->_ndimension - coarse->_ndimension;
+    assert(lowerdims >= 0);
+    for(int d=0;d<lowerdims;d++){
+      assert(fine->_simd_layout[d]==1);
+      assert(fine->_processors[d]==1);
+    }
 
-    multiplicity = multiplicity *fine->_rdimensions[fd] / coarse->_rdimensions[d]; 
+    int multiplicity=1;
+    for(int d=0;d<lowerdims;d++){
+      multiplicity=multiplicity*fine->_rdimensions[d];
+    }
+    // local and global volumes subdivide cleanly after SIMDization
+    for(int d=0;d<rngdims;d++){
+      int fd= d+lowerdims;
+      assert(coarse->_processors[d]  == fine->_processors[fd]);
+      assert(coarse->_simd_layout[d] == fine->_simd_layout[fd]);
+      assert(((fine->_rdimensions[fd] / coarse->_rdimensions[d])* coarse->_rdimensions[d])==fine->_rdimensions[fd]); 
+
+      multiplicity = multiplicity *fine->_rdimensions[fd] / coarse->_rdimensions[d]; 
+    }
+    return multiplicity;
   }
-  return multiplicity;
 }
 
   
@@ -80,6 +94,19 @@ inline int RNGfillable(GridBase *coarse,GridBase *fine)
 // this function is necessary for the LS vectorised field
 inline int RNGfillable_general(GridBase *coarse,GridBase *fine)
 {
+
+  if ( coarse == fine ) return 1;
+
+  if ( coarse->isIcosahedral()) assert(coarse->isIcosahedralEdge());
+  
+  if ( fine->isIcosahedralVertex() && coarse->isIcosahedralEdge() ) {
+    assert(fine->Nd()==coarse->Nd());
+    for(int d=0;d<fine->Nd();d++){
+      assert(fine->LocalDimensions()[d] == coarse->LocalDimensions()[d]);
+    }
+    return 1;
+  }
+
   int rngdims = coarse->_ndimension;
     
   // trivially extended in higher dims, with locality guaranteeing RNG state is local to node
@@ -352,12 +379,12 @@ private:
 public:
   GridBase *Grid(void) const { return _grid; }
   int generator_idx(int os,int is) {
-    return is*_grid->oSites()+os;
+    return (is*_grid->CartesianOsites()+os)%_grid->lSites(); // On the pole sites wrap back to normal generators; Icosahedral hack
   }
 
   GridParallelRNG(GridBase *grid) : GridRNGbase() {
     _grid = grid;
-    _vol  =_grid->iSites()*_grid->oSites();
+    _vol  =_grid->lSites();
 
     _generators.resize(_vol);
     _uniform.resize(_vol,std::uniform_real_distribution<RealD>{0,1});
@@ -381,7 +408,7 @@ public:
 
     int multiplicity = RNGfillable_general(_grid, l.Grid()); // l has finer or same grid
     int Nsimd  = _grid->Nsimd();  // guaranteed to be the same for l.Grid() too
-    int osites = _grid->oSites();  // guaranteed to be <= l.Grid()->oSites() by a factor multiplicity
+    int osites = _grid->CartesianOsites();  // guaranteed to be <= l.Grid()->oSites() by a factor multiplicity, except on Icosahedral
     int words  = sizeof(scalar_object) / sizeof(scalar_type);
 
     autoView(l_v, l, CpuWrite);
@@ -402,8 +429,27 @@ public:
 	// merge into SIMD lanes, FIXME suboptimal implementation
 	merge(l_v[sm], buf);
       }
-      });
-    //    });
+    });
+
+    /*
+     * Fill in the poles for an Icosahedral vertex mesh
+     */
+    if (l.Grid()->isIcosahedralVertex()) { 
+      int64_t pole_sites=l.Grid()->NorthPoleOsites()+l.Grid()->SouthPoleOsites();
+      int64_t pole_base =l.Grid()->CartesianOsites();
+
+      ExtractBuffer<scalar_object> buf(Nsimd);
+      for (int m = 0; m < pole_sites; m++) {  // Draw from same generator multiplicity times                                                                                                           
+        for (int si = 0; si < Nsimd; si++) {
+          int gdx = 0;
+	  scalar_type *pointer = (scalar_type *)&buf[si];
+          dist[gdx].reset();
+          for (int idx = 0; idx < words; idx++)
+            fillScalar(pointer[idx], dist[gdx], _generators[gdx]);
+        }
+        merge(l_v[pole_base+m], buf);
+      }      
+    }
 
     _time_counter += usecond()- inner_time_counter;
   }

Grid/stencil/IcosahedralStencil.h

@@ -0,0 +1,768 @@
+/*************************************************************************************
+
+     Grid physics library, www.github.com/paboyle/Grid 
+
+     Source file: ./lib/GeneralLocalStencil.h
+
+     Copyright (C) 2019
+
+ Author: Peter Boyle <paboyle@ph.ed.ac.uk>
+
+     This program is free software; you can redistribute it and/or modify
+     it under the terms of the GNU General Public License as published by
+     the Free Software Foundation; either version 2 of the License, or
+     (at your option) any later version.
+
+     This program is distributed in the hope that it will be useful,
+     but WITHOUT ANY WARRANTY; without even the implied warranty of
+     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+     GNU General Public License for more details.
+
+     You should have received a copy of the GNU General Public License along
+     with this program; if not, write to the Free Software Foundation, Inc.,
+     51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+
+     See the full license in the file "LICENSE" in the top level distribution directory
+*************************************************************************************/
+/*  END LEGAL */
+#pragma once
+NAMESPACE_BEGIN(Grid);
+
+// Share with Cartesian Stencil
+struct IcosahedralStencilEntry { 
+  uint64_t _offset;            // 8 bytes 
+  uint8_t _is_local;           // is this with our lattice array (else in a comms buf)
+  uint8_t _adjoint;            // is this with our lattice array (else in a comms buf)
+  uint8_t _polarisation;       // which lorentz index on the neighbours patch
+  uint8_t _missing_link;       // 
+};
+
+enum IcoshedralDirections {
+  IcosahedronPatchX = 0,
+  IcosahedronPatchY = 1,
+  IcosahedronPatchDiagonal=2,
+  IcosahedronTime=3
+};
+
+inline int periAdd(int A,int inc,int L) { return (A+inc+L)%L ; }
+
+class IcosahedralStencilView {
+ public:
+  ////////////////////////////////////////
+  // Basic Grid and stencil info
+  ////////////////////////////////////////
+  int                       _npoints; // Move to template param?
+  IcosahedralStencilEntry*  _entries_p;
+
+  accelerator_inline IcosahedralStencilEntry * GetEntry(int point,int osite) const { 
+    return & this->_entries_p[point+this->_npoints*osite]; 
+  }
+  void ViewClose(void){};
+};
+
+////////////////////////////////////////
+// The Stencil Class itself
+////////////////////////////////////////
+class IcosahedralStencil : public IcosahedralStencilView {
+public:
+  typedef IcosahedralStencilView View_type;
+
+protected:
+  GridBase *                        _grid;
+
+public: 
+  GridBase *Grid(void) const { return _grid; }
+
+  View_type View(int mode) const {
+    View_type accessor(*( (View_type *) this));
+    return accessor;
+  }
+  // NB x+, y+ is ALWAYS present, as are the forward 3 directions for links owned by each site
+  //
+  // These are VERTEX mesh neigbours being returned, with isPole indicating we need N/S pole according to
+  // hemisphere.
+  //
+  // If needing edge mesh "neigbours" to assemble loops we must find the mapping of a forward link
+  // to a corresponding "backward link" on the pole
+  deviceVector<IcosahedralStencilEntry>  _entries;
+  
+  void GetNbrForPlusX(GridBase *grid,Coordinate &Coor,Coordinate &NbrCoor, int &isPole)
+  {
+    int nd = grid->Nd();
+    int L  = grid->LocalDimensions()[0];
+    assert (grid->LocalDimensions()[0] == grid->LocalDimensions()[1]);
+    int patch = Coor[nd-1];
+    int north = patch/HemiPatches;
+    int south = 1-north;
+    int HemiPatch = patch%HemiPatches;
+    int NbrPatch;
+    
+    NbrCoor = Coor;
+    isPole = 0;
+    if ( Coor[0]<(L-1) ) { 
+      NbrCoor[0]=Coor[0]+1;
+      // Nbr is inside THIS patch
+      return;
+    }
+    if ( north ) {
+      assert(Coor[0]==(L-1));
+      // Simple connect to southern neighbour tile
+      NbrCoor[0]=0;
+      NbrCoor[nd-1]=periAdd(HemiPatch,+1,HemiPatches) + SouthernHemisphere;
+      return;
+    }
+    ////////////////////////////////////////////////////////////
+    // FIXME:
+    // Can store the rotation of polarisation here: get xdir instead of ydir and must take adjoint
+    ////////////////////////////////////////////////////////////
+    if ( south ) {
+      assert(Coor[0]==(L-1));
+      if ( Coor[1] == 0 ) {
+	isPole = 1;
+	NbrCoor[0] = (L-1); // Coordinate of the edge graph site holding the edge for other vertex in pole triangle
+	NbrCoor[1] = 0;
+	NbrCoor[nd-1]=periAdd(HemiPatch,+1,HemiPatches) + SouthernHemisphere;
+	return;
+      } else {
+	NbrCoor[1] = 0;
+	NbrCoor[0] = L-Coor[1];
+	NbrCoor[nd-1]=periAdd(HemiPatch,+1,HemiPatches) + SouthernHemisphere;
+	return;
+      }
+    }
+    assert(0);
+  }
+  void GetNbrForPlusY(GridBase *grid,Coordinate &Coor,Coordinate &NbrCoor, int &isPole)
+  {
+    int nd = grid->Nd();
+    int L  = grid->LocalDimensions()[0];
+    assert (grid->LocalDimensions()[0] == grid->LocalDimensions()[1]);
+    int patch = Coor[nd-1];
+    int north = patch/HemiPatches;
+    int south = 1-north;
+    int HemiPatch = patch%HemiPatches;
+    int NbrPatch;
+    
+    NbrCoor = Coor;
+    isPole = 0;
+    if ( Coor[1]<(L-1) ) { 
+      NbrCoor[1]=Coor[1]+1;
+      // Nbr is inside THIS patch
+      return;
+    }
+    if ( south ) {
+      // Simple connect to northern neighbour tile
+      assert(Coor[1]==(L-1));
+      NbrCoor[1]=0;
+      NbrCoor[nd-1]=HemiPatch + NorthernHemisphere;
+      return;
+    }
+    ////////////////////////////////////////////////////////////
+    // FIXME:
+    // Could store the rotation of polarisation here: get xdir instead of ydir and must take adjoint
+    // But probaby just write "getLinkPropertiesToCloseTriangleA" "getLinkPropertiesToCloseTriangleB"
+    // Or write a 'double store' method to move edge graph to a vertex graph with all fermion transports
+    ////////////////////////////////////////////////////////////
+    if ( north ) {
+      assert(Coor[1]==(L-1));
+      if ( Coor[0] == 0 ) {
+	isPole = 1;
+	NbrCoor[1] = (L-1); // Coordinate of the edge graph site holding the edge for other vertex in pole triangle
+	NbrCoor[0] = 0;
+	NbrCoor[nd-1]=periAdd(HemiPatch,+1,HemiPatches) + NorthernHemisphere;
+	return;
+      } else {
+	NbrCoor[0] = 0;
+	NbrCoor[1] = L-Coor[0]; // x=1 --> y=L-1 for y+
+	NbrCoor[nd-1]=periAdd(HemiPatch,+1,HemiPatches) + NorthernHemisphere;
+	return;
+      }
+    }
+    assert(0);
+  }
+  // Missing links are at (0,0) on local patch coordinates in -diagonal direction 
+  // We are here returning VERTEX grid coordinates.
+  void GetNbrForMinusX(GridBase *grid,Coordinate &Coor,Coordinate &NbrCoor)
+  {
+    int nd = grid->Nd();
+    int L  = grid->LocalDimensions()[0];
+    assert (grid->LocalDimensions()[0] == grid->LocalDimensions()[1]);
+    int patch = Coor[nd-1];
+    int north = patch/HemiPatches;
+    int south = 1-north;
+    int HemiPatch = patch%HemiPatches;
+    int NbrPatch;
+    
+    NbrCoor = Coor;
+    if ( Coor[0]>0 ) { 
+      NbrCoor[0]=Coor[0]-1;
+      return;
+    }
+    if ( south ) {
+      assert(Coor[0]==0);
+      // Simple connect to northern neighbour tile
+      NbrCoor[0]=L-1;
+      NbrCoor[nd-1]=periAdd(HemiPatch,-1,HemiPatches) + NorthernHemisphere;
+      return;
+    }
+    if ( north ) {
+	NbrCoor[0] = L-1-Coor[1];
+	NbrCoor[1] = L-1;
+	NbrCoor[nd-1]=periAdd(HemiPatch,+1,HemiPatches) + SouthernHemisphere;
+	return;
+    }      
+    assert(0);
+  }
+  
+  void GetNbrForMinusY(GridBase *grid,Coordinate &Coor,Coordinate &NbrCoor)
+  {
+    int nd = grid->Nd();
+    int L  = grid->LocalDimensions()[0];
+    assert (grid->LocalDimensions()[0] == grid->LocalDimensions()[1]);
+    int patch = Coor[nd-1];
+    int north = patch/HemiPatches;
+    int south = 1-north;
+    int HemiPatch = patch%HemiPatches;
+    int NbrPatch;
+    
+    NbrCoor = Coor;
+    if ( Coor[1]>0 ) { 
+      NbrCoor[1]=Coor[1]-1;
+      return;
+    }
+    if ( north ) {
+      assert(Coor[1]==0);
+      // Simple connect to northern neighbour tile
+      NbrCoor[1]=L-1;
+      NbrCoor[nd-1]=HemiPatch + SouthernHemisphere;
+      return;
+    }
+    if ( south ) {
+	NbrCoor[1] = L-1-Coor[0];
+	NbrCoor[0] = L-1;
+	NbrCoor[nd-1]=periAdd(HemiPatch,-1,HemiPatches) + SouthernHemisphere;
+	return;
+    }      
+    assert(0);
+  }
+  void GetNbrForMinusDiagonal(GridBase *grid,Coordinate &Coor,Coordinate &NbrCoor,int &missingLink)
+  {
+    int nd = grid->Nd();
+    int L  = grid->LocalDimensions()[0];
+    assert (grid->LocalDimensions()[0] == grid->LocalDimensions()[1]);
+    int patch = Coor[nd-1];
+    int north = patch/HemiPatches;
+    int south = 1-north;
+    int HemiPatch = patch%HemiPatches;
+    int NbrPatch;
+    
+    NbrCoor = Coor;
+    
+    missingLink = 0;
+    if( Coor[0]==0 && Coor[1]==0) {
+      missingLink=1;
+      return;
+    }
+    if ( (Coor[0]>0) && (Coor[1]>0) ) { 
+      // Nbr is inside THIS patch
+      NbrCoor[0]=Coor[0]-1;
+      NbrCoor[1]=Coor[1]-1;
+      return;
+    }
+    if ( north ) {
+      if ( Coor[0]==0 ) {
+	// We are on -x edge
+	// Maps to +y edge of hemipatch to the left
+	NbrCoor[nd-1] = periAdd(HemiPatch,-1,HemiPatches) + NorthernHemisphere;
+	NbrCoor[0]=(L-Coor[1]);
+	NbrCoor[1]=(L-1);
+	return;
+      } else {
+        // We are on the -y edge and NOT bottom left corner; Nbr is in the patch LEFT
+	assert( (Coor[0]>0) && (Coor[1]==0) );
+	NbrCoor[nd-1] = HemiPatch + SouthernHemisphere; // Map from north to south
+	NbrCoor[0]=Coor[0]-1;
+	NbrCoor[1]=(L-1);
+	return;
+      }
+      assert(0);
+    }
+    if ( south ) {
+      // We are on the -y edge
+      if ( Coor[1]==0 ) {
+	NbrCoor[nd-1] = periAdd(HemiPatch,-1,HemiPatches) + SouthernHemisphere;
+        NbrCoor[0]=(L-1);
+        NbrCoor[1]=(L-Coor[0]);
+	return;
+      } else {
+      // We are on the -x edge and NOT bottom left corner; Nbr is in the patch LEFT
+	assert( (Coor[0]==0) && (Coor[1]>0) );
+	NbrCoor[nd-1] = periAdd(HemiPatch,-1,HemiPatches) + NorthernHemisphere; // south to north
+	NbrCoor[0]=(L-1);
+	NbrCoor[1]= Coor[1]-1;
+	return;
+      }
+      assert(0);
+    }
+    assert(0);
+  }
+  void GetNbrForPlusDiagonal(GridBase *grid,Coordinate &Coor,Coordinate &NbrCoor)
+  {
+    int nd = grid->Nd();
+    int L  = grid->LocalDimensions()[0];
+    assert (grid->LocalDimensions()[0] == grid->LocalDimensions()[1]);
+    int patch = Coor[nd-1];
+    int north = patch/HemiPatches;
+    int south = 1-north;
+    int HemiPatch = patch%HemiPatches;
+    int NbrPatch;
+    
+    NbrCoor = Coor;
+    
+    if ( (Coor[0]<L-1) && (Coor[1]<L-1) ) { 
+      // Nbr is inside THIS patch
+      NbrCoor[0]=Coor[0]+1; 
+      NbrCoor[1]=Coor[1]+1;
+      return;
+    }
+
+    if ( north ) {
+      // We are on +y edge
+      if ( Coor[1]==(L-1) ) {
+	NbrCoor[nd-1] = periAdd(HemiPatch,+1,HemiPatches) + NorthernHemisphere;
+	NbrCoor[0]=0;
+	NbrCoor[1]=(L-1)-Coor[0];
+	return;
+      } else {
+      // Else we are on the +x edge, not top right corner
+	assert( Coor[0]==(L-1) && Coor[1]<(L-1) );
+	NbrCoor[nd-1] = periAdd(HemiPatch,+1,HemiPatches) + SouthernHemisphere;
+	NbrCoor[0]=0;
+	NbrCoor[1]=Coor[1]+1;
+	return;
+      }
+      assert(0);
+    }
+
+    if ( south ) {
+      // We are on the +x edge
+      if ( Coor[0]==(L-1) ) {
+	NbrCoor[nd-1] = periAdd(HemiPatch,+1,HemiPatches) + SouthernHemisphere;
+	NbrCoor[0]=(L-1)-Coor[1]; //y=(L-1) <-> x=0 ; y=0<->x=(L-1)  [ Sanity check ]
+	NbrCoor[1]=0;
+	return;
+      } else { 
+      // We are on the +y edge and NOT top right corner; Nbr is in the patch UP
+	assert( (Coor[1]==L-1) && (Coor[0]<(L-1)) );
+	NbrCoor[nd-1] = HemiPatch + NorthernHemisphere;
+	NbrCoor[0]=Coor[0]+1;
+	NbrCoor[1]=0;
+	return;
+      }
+      assert(0);
+    }
+    assert(0);
+  }
+  
+  void  TestGeometry(void)
+  {
+    GridBase *grid = this->_grid;
+    uint64_t cart_sites = grid->CartesianOsites();
+
+    //////////////////////////////////////////////////////////////////////////////////////////////
+    //    Loop over cart sites.
+    //    Find two triangles per site.
+    //    Check going forward in X, Up and forward in Diag match
+    //    Check going Up, forward in X and forward Diag match; subtleties at poles and rotation in cross patch
+    //////////////////////////////////////////////////////////////////////////////////////////////
+
+    std::cout << "*************************************"<<std::endl;
+    std::cout << " Icosahedral Stencil Geometry Test !"<<std::endl;
+    std::cout << "*************************************"<<std::endl;
+
+    const int triangle_ref = cart_sites;
+    std::cout << " Base triangle count for each type " <<triangle_ref;
+    
+    std::cout << "------------------------------------"<<std::endl;
+    std::cout << "testing +x+y vs +diag"<<std::endl;
+    std::cout << "testing +y+x vs +diag"<<std::endl;
+    std::cout << "------------------------------------"<<std::endl;
+    int xyd_pole_count=0;
+    int xyd_count=0;
+    int yxd_pole_count=0;
+    int yxd_count=0;
+    for(uint64_t site=0;site<cart_sites; site ++) {
+
+      Coordinate Coor;
+      Coordinate DiagCoor;
+
+      int nd = grid->Nd();
+      int L  = grid->LocalDimensions()[0];
+
+      grid->oCoorFromOindex(Coor,site);
+
+      int patch = Coor[nd-1];
+      int north = patch/HemiPatches;
+      int south = 1-north;
+      int isPole      = 0;
+      int discard;
+      int missingLink = 0;
+      
+      int HemiPatch = patch%HemiPatches;
+
+      //////////////////////////////
+      // First test of triangle
+      //////////////////////////////
+      // Compare +x, +y to +diag
+      Coordinate XpCoor;
+      Coordinate YpXpCoor;
+      GetNbrForPlusDiagonal(grid,Coor,DiagCoor);
+      GetNbrForPlusX(grid,Coor,XpCoor,isPole);
+
+      int XpHemiPatch  = XpCoor[nd-1]%HemiPatches;
+      int XpHemisphere = XpCoor[nd-1]/HemiPatches;
+      
+      if (isPole) {
+	YpXpCoor = XpCoor;
+      } else if ( XpHemiPatch != HemiPatch && south ) {
+	GetNbrForMinusX(grid,XpCoor,YpXpCoor);
+      } else {
+	GetNbrForPlusY(grid,XpCoor,YpXpCoor,discard);
+      }
+      
+      if(isPole) {
+	std::cout << "Forward xyd triangle "<<Coor<<"-Pole["<<XpCoor[2]<<"]-"<<YpXpCoor<<" should be " <<DiagCoor<<std::endl;
+	xyd_pole_count++;
+      } else {
+	std::cout << "Forward xyd triangle "<<Coor<<"-"<<XpCoor<<"-"<<YpXpCoor<<" should be " <<DiagCoor<<std::endl;
+	xyd_count++;
+      }
+      for(int d=0;d<DiagCoor.size();d++) {
+	assert(DiagCoor[d]==YpXpCoor[d]);
+      }
+
+      Coordinate YpCoor;
+      Coordinate XpYpCoor;
+      GetNbrForPlusDiagonal(grid,Coor,DiagCoor);
+      GetNbrForPlusY(grid,Coor,YpCoor,isPole);
+
+      int YpHemiPatch  = YpCoor[nd-1]%HemiPatches;
+      int YpHemisphere = YpCoor[nd-1]/HemiPatches;
+      
+      if(isPole) {
+	XpYpCoor = YpCoor;
+      } else if ( YpHemiPatch != HemiPatch  && north ) {
+	GetNbrForMinusY(grid,YpCoor,XpYpCoor); // we hopped - this rotates the directions
+      } else {
+	GetNbrForPlusX(grid,YpCoor,XpYpCoor,discard);
+      }
+
+      if(isPole) {
+	yxd_pole_count++;
+	std::cout << "Forward yxd triangle "<<Coor<<"-Pole["<<YpCoor[2]<<"]-"<<XpYpCoor<<" should be " <<DiagCoor<<std::endl;
+      } else {
+	yxd_count++;
+	std::cout << "Forward yxd triangle "<<Coor<<"-"<<YpCoor<<"-"<<XpYpCoor<<" should be " <<DiagCoor<<std::endl;
+      }
+      for(int d=0;d<DiagCoor.size();d++) {
+	assert(DiagCoor[d]==XpYpCoor[d]);
+      }
+    }
+    std::cout << " xyd_count "<<xyd_count<<" + poles_count "<<xyd_pole_count<<" expect "<<triangle_ref<<" triangles "<<std::endl;
+    std::cout << " yxd_count "<<yxd_count<<" + poles_count "<<yxd_pole_count<<" expect "<<triangle_ref<<" triangles "<<std::endl;
+    assert(xyd_count+xyd_pole_count == triangle_ref);
+    assert(yxd_count+yxd_pole_count == triangle_ref);
+    std::cout << "------------------------------------"<<std::endl;
+    std::cout << "testing -diag +x+y = identity"<<std::endl;
+    std::cout << "testing -diag +y+x = identity"<<std::endl;
+    std::cout << "------------------------------------"<<std::endl;
+
+    int dmxy_count=0;
+    int dmyx_count=0;
+    int dmxy_count_special=0;
+    int dmyx_count_special=0;
+    int num_missing=0;
+
+    for(uint64_t site=0;site<cart_sites; site ++) {
+
+      Coordinate Coor;
+
+      int nd = grid->Nd();
+      int L  = grid->LocalDimensions()[0];
+
+      grid->oCoorFromOindex(Coor,site);
+
+      int patch = Coor[nd-1];
+      int north = patch/HemiPatches;
+      int south = 1-north;
+      int isPole      = 0;
+      int discard;
+      int missingLink = 0;
+      int HemiPatch = patch%HemiPatches;
+
+      Coordinate DmCoor;
+      GetNbrForMinusDiagonal(grid,Coor,DmCoor,missingLink);
+      if ( missingLink ) {
+	std::cout << Coor << " has no backwards diagonal link "<<std::endl;
+	num_missing++;
+      } else {
+
+	int DmPatch  = DmCoor[nd-1];
+	int DmHemiPatch  = DmCoor[nd-1]%HemiPatches;
+	int DmHemisphere = DmCoor[nd-1]/HemiPatches;
+
+	Coordinate XpDmCoor;
+	Coordinate YpXpDmCoor;
+
+	Coordinate YpDmCoor;
+	Coordinate XpYpDmCoor;
+	
+	if ( DmHemiPatch != HemiPatch && north ) {
+	  GetNbrForPlusDiagonal(grid,DmCoor,XpDmCoor);
+	  GetNbrForPlusY(grid,XpDmCoor,YpXpDmCoor,isPole); assert(!isPole);
+
+	  GetNbrForMinusX(grid,DmCoor,YpDmCoor);
+	  GetNbrForPlusDiagonal(grid,YpDmCoor,XpYpDmCoor);
+	  dmxy_count_special++;
+	  dmyx_count_special++;
+	} else if ( DmPatch == periAdd(HemiPatch,-1,HemiPatches) && south ) {
+	  GetNbrForPlusDiagonal(grid,DmCoor,YpDmCoor);
+	  GetNbrForPlusX(grid,YpDmCoor,XpYpDmCoor,isPole); assert(!isPole);
+
+	  GetNbrForMinusY(grid,DmCoor,XpDmCoor);
+	  GetNbrForPlusDiagonal(grid,XpDmCoor,YpXpDmCoor);
+	  dmxy_count_special++;
+	  dmyx_count_special++;
+	} else {
+	  GetNbrForPlusX(grid,DmCoor,XpDmCoor,isPole); assert(!isPole);
+	  GetNbrForPlusY(grid,XpDmCoor,YpXpDmCoor,isPole); assert(!isPole);
+
+	  GetNbrForPlusY(grid,DmCoor,YpDmCoor,isPole); assert(!isPole);
+	  GetNbrForPlusX(grid,YpDmCoor,XpYpDmCoor,isPole); assert(!isPole);
+	  dmxy_count++;
+	  dmyx_count++;
+	}
+	std::cout << Coor<<" DmXpYp triangle YpXpDm"<<YpXpDmCoor<<"-XpDm"<<XpDmCoor<<"-Dm"<<DmCoor<<" should be " <<Coor<<std::endl;
+	for(int d=0;d<Coor.size();d++) {
+	  assert(Coor[d]==YpXpDmCoor[d]);
+	}
+
+	std::cout << Coor<<"DmXpYp triangle XpYpDm"<<XpYpDmCoor<<"-YpDm"<<YpDmCoor<<"-Dm"<<DmCoor<<" should be " <<Coor<<std::endl;
+	for(int d=0;d<Coor.size();d++) {
+	  assert(Coor[d]==XpYpDmCoor[d]);
+	}
+      }
+    }
+    std::cout << " dmxy_count "<<dmxy_count<<" + special "<<dmxy_count_special<<" + missing "<<num_missing<<" expect "<<triangle_ref<<" triangles "<<std::endl;
+      std::cout << " dmyx_count "<<dmyx_count<<" + special "<<dmyx_count_special<<" + missing "<<num_missing<<" expect "<<triangle_ref<<" triangles "<<std::endl;
+    assert(dmxy_count + dmxy_count_special + num_missing == triangle_ref);
+    assert(dmyx_count + dmyx_count_special + num_missing == triangle_ref);
+    
+    std::cout << "------------------------------------"<<std::endl;
+    std::cout << "testing diag -x-y = identity "<<std::endl;
+    std::cout << "testing diag -y-x = identity"<<std::endl;
+    std::cout << "------------------------------------"<<std::endl;
+
+    std::cout << "------------------------------------"<<std::endl;
+    std::cout << "testing -diag = -x-y "<<std::endl;
+    std::cout << "testing -diag = -y-x "<<std::endl;
+    std::cout << "------------------------------------"<<std::endl;
+    
+    std::cout << "*************************************"<<std::endl;
+    std::cout << " Icosahedral Stencil Geometry Test Complete"<<std::endl;
+    std::cout << "*************************************"<<std::endl;
+  }
+  IcosahedralStencil(GridBase *grid) // Must be +1 or -1
+  {
+    this->_grid    = grid;
+    // Loop over L^2 x T x npatch and the
+    assert(grid->isIcosahedral());
+  }
+  void NearestNeighbourStencil(void)
+  {
+    GridBase * grid = this->_grid;
+    assert(grid->isIcosahedralVertex());
+    
+  }
+
+  
+  /*
+   * For gauge action implementation
+   */
+  void FaceStencil(void)
+  {
+    GridBase * grid = this->_grid;
+
+    int osites  = grid->oSites();
+
+    uint64_t cart_sites = grid->CartesianOsites();
+    uint64_t Npole_sites = grid->NorthPoleOsites();
+    uint64_t Spole_sites = grid->SouthPoleOsites();
+    Coordinate pcoor = grid->ThisProcessorCoor();
+    Coordinate pgrid = grid->ProcessorGrid();
+
+    /*
+     * resize the stencil entries array and set npoints
+     */
+    this->_npoints=2; // Move to template param?
+    this->_entries.resize(this->_npoints * cart_sites);
+    this->_entries_p = &_entries[0];
+
+    for(uint64_t site=0;site<cart_sites; site ++) {
+
+      Coordinate Coor;
+      Coordinate NbrCoor;
+
+      int nd = grid->Nd();
+      int L  = grid->LocalDimensions()[0];
+
+      Integer lexXY = site*2;
+      Integer lexYX = site*2+1;
+
+      IcosahedralStencilEntry SE;
+	
+      ////////////////////////////////////////////////
+      // Outer index of neighbour Offset calculation
+      ////////////////////////////////////////////////
+      grid->oCoorFromOindex(Coor,site);
+      NbrCoor = Coor;
+      assert( grid->LocalDimensions()[1]==grid->LocalDimensions()[0]);
+      assert( grid->_simd_layout[0]==1); // Cannot vectorise in these dims
+      assert( grid->_simd_layout[1]==1);
+      assert( grid->_processors[0]==1); // Cannot mpi distribute in these dims
+      assert( grid->_processors[1]==1);
+
+      int Patch = Coor[nd-1];
+      int HemiPatch = Patch%HemiPatches;
+      int north = Patch/HemiPatches;
+      int south = 1-north;
+      int isPoleY;
+      int isPoleX;
+      
+      assert(Patch<IcosahedralPatches);
+      assert((north==1)||(south==1));
+	/*
+	 * Minimal stencil for edge -> face triangle evaluation
+	 *
+	 * On each edge grid site, hold +x,+y,+diag links
+	 * Must locate the closing link to form the two forward triangles
+	 *
+	 * case: +x,+d 
+	 *
+	 *   north: +x neighbours ydir link always
+	 *   south: +x neighbours ydir link except
+	 *            cross into a different patch in south
+	 *         then
+	 *          +d neighbours xdir link
+	 *
+	 * case: +y,+d 
+	 *   south: +y neighbours xdir link always
+	 *   north: +y neighbours xdir link unless
+	 *            cross into a different patch in north
+	 *         then
+	 *          +d neighbours ydir link
+	 *
+	 */
+      Coordinate XpCoor;
+      Coordinate YpCoor;
+      Coordinate DpCoor;
+
+      GetNbrForPlusDiagonal(grid,Coor,DpCoor);
+      GetNbrForPlusX(grid,Coor,XpCoor,isPoleX);
+      GetNbrForPlusY(grid,Coor,YpCoor,isPoleY);
+
+      int XpHemiPatch  = XpCoor[nd-1]%HemiPatches;
+      int XpHemisphere = XpCoor[nd-1]/HemiPatches;
+
+      int DpPatch  = DpCoor[nd-1];
+      int DpHemiPatch  = DpCoor[nd-1]%HemiPatches;
+      int DpHemisphere = DpCoor[nd-1]/HemiPatches;
+
+      ////////////////////////////////////////////////
+      // for trace [ U_x(z) U_y(z+\hat x) adj(U_d(z)) ]
+      ////////////////////////////////////////////////
+      if ( DpHemiPatch != HemiPatch && south ) {
+	SE._offset       = grid->oIndex(DpCoor);
+	SE._is_local     = true;
+	SE._polarisation = IcosahedronPatchX;
+	SE._adjoint      = true;
+	SE._missing_link = false;
+      } else {
+	SE._offset       = grid->oIndex(XpCoor);
+	SE._is_local     = true;
+	SE._polarisation = IcosahedronPatchY;
+	SE._adjoint      = false;
+	SE._missing_link = false;
+      }
+
+      ////////////////////////////////////////////////
+      // Store in look up table
+      ////////////////////////////////////////////////
+      acceleratorPut(this->_entries[lexXY],SE);
+      
+      ////////////////////////////////////////////////
+      // for trace [  U_y(z) adj(U_d(z))  U_x(z+\hat y) ]
+      ////////////////////////////////////////////////
+      if ( DpHemiPatch != HemiPatch && north ) {
+	SE._offset       = grid->oIndex(DpCoor);
+	SE._is_local     = true;
+	SE._polarisation = IcosahedronPatchY;
+	SE._adjoint      = true;
+	SE._missing_link = false;
+      } else {
+	SE._offset       = grid->oIndex(YpCoor);
+	SE._is_local     = true;
+	SE._polarisation = IcosahedronPatchX;
+	SE._adjoint      = false;
+	SE._missing_link = false;
+      }
+      ////////////////////////////////////////////////
+      // Store in look up table
+      ////////////////////////////////////////////////
+      acceleratorPut(this->_entries[lexYX],SE);
+    };
+  }
+
+  /*
+   * For gauge action derivative implementation
+   * Staple 
+   *
+   * Case1: I x T loops
+   *
+   * Need: DirP this site; no entry
+   *       Tp   @ dir++
+   *       DirP @ t++
+   *       Tp   @ t--
+   *       DirP @ t--
+   *       Tp   @ t--, dir++
+   *
+   * There is no complex rotation of links on other site
+   *
+   * Case2: I x I loops
+   *
+   * Y staple: need
+   *     Diag @ (xy)
+   *     X    @ y++  ; care needed for rotation
+   *     Diag @ x--  ; care needed for rotation
+   *     X    @ x--  ; care needed for rotation
+   *
+   * X staple: need
+   *     Diag @ (xy)
+   *     X    @ y++  ; care needed for rotation
+   *     Diag @ x--  ; care needed for rotation
+   *     X    @ x--  ; care needed for rotation
+   *
+   * Diag staple: need
+   *
+   *     X@  (xy)
+   *     Y@  x++  ; care needed for rotation
+   *     Y@  (xy)
+   *     X@  y++  ; care needed for rotation
+   */
+  void StapleDiagStencil(void){  }
+  void StapleXpStencil(void)  {  }
+  void StapleYpStencil(void)  {  }
+  void StapleTpStencil(void)  {  }
+  
+};
+NAMESPACE_END(Grid);

systems/mac-arm/config-command

@@ -1,5 +1,6 @@
 MPICXX=mpicxx CXXFLAGS=-I/opt/local/include LDFLAGS=-L/opt/local/lib/ CXX=clang++ ../../configure \
 	      --enable-simd=GEN \
+	      --enable-Nc=1 \
 	      --enable-comms=mpi-auto \
 	      --enable-unified=yes \
 	      --prefix $HOME/QCD/GridInstall \

tests/debug/Test_icosahedron.cc

@@ -26,27 +26,231 @@ Author: Peter Boyle <paboyle@ph.ed.ac.uk>
     *************************************************************************************/
     /*  END LEGAL */
 #include <Grid/Grid.h>
+#include <Grid/stencil/IcosahedralStencil.h>
 
 using namespace std;
 using namespace Grid;
 
 const int MyNd=3;
-template<typename vtype> using iIcosahedralLorentzComplex = iVector<iScalar<iScalar<vtype> >, MyNd+1 > ;
 
-typedef iIcosahedralLorentzComplex<ComplexD >  IcosahedralLorentzComplexD;
-typedef iIcosahedralLorentzComplex<vComplexD> vIcosahedralLorentzComplexD;
-typedef Lattice<vIcosahedralLorentzComplexD> LatticeIcosahedralLorentzComplexD;
+template<typename vtype> using iIcosahedralLorentzComplex      = iVector<iScalar<iScalar<vtype> >, MyNd+1 > ;
+template<typename vtype> using iIcosahedralLorentzColourMatrix = iVector<iScalar<iMatrix<vtype,Nc> >, MyNd+1 > ;
+template<typename vtype> using iIcosahedralColourMatrix        = iScalar<iScalar<iMatrix<vtype,Nc> > > ;
+
+typedef iIcosahedralLorentzComplex<Complex >  IcosahedralLorentzComplex;
+typedef iIcosahedralLorentzComplex<vComplex> vIcosahedralLorentzComplex;
+typedef Lattice<vIcosahedralLorentzComplex> LatticeIcosahedralLorentzComplex;
+
+typedef iIcosahedralLorentzColourMatrix<Complex >  IcosahedralLorentzColourMatrix;
+typedef iIcosahedralLorentzColourMatrix<vComplex> vIcosahedralLorentzColourMatrix;
+typedef Lattice<vIcosahedralLorentzColourMatrix>  LatticeIcosahedralLorentzColourMatrix;
+
+typedef iIcosahedralColourMatrix<Complex >  IcosahedralColourMatrix;
+typedef iIcosahedralColourMatrix<vComplex> vIcosahedralColourMatrix;
+typedef Lattice<vIcosahedralColourMatrix>  LatticeIcosahedralColourMatrix;
+
+
+class IcosahedralGimpl
+{
+public:
+  typedef LatticeIcosahedralLorentzColourMatrix  GaugeField;
+  typedef LatticeIcosahedralColourMatrix         GaugeLinkField;
+  typedef LatticeComplex                          ComplexField;
+};
+  
+template< class Gimpl>
+class IcosahedralEdgeSupport {
+public:
+  // Move to inherit types macro as before
+  typedef typename Gimpl::GaugeField GaugeField;
+  typedef typename Gimpl::GaugeLinkField GaugeLinkField;
+  typedef typename Gimpl::ComplexField ComplexField;
+  //
+  GridBase *VertexGrid;
+  GridBase *EdgeGrid;
+  IcosahedralStencil FaceStencil;
+  
+  IcosahedralEdgeSupport(GridBase *_VertexGrid,GridBase *_EdgeGrid)
+    : FaceStencil (EdgeGrid), VertexGrid(_VertexGrid), EdgeGrid(_EdgeGrid)
+  {
+    FaceStencil.FaceStencil();
+  }
+
+  ////////////////////////////////////////////////////////////////////////////////////
+  // Fixme: will need a version of "Gimpl" and a wrapper class following "WilsonLoops" style.
+  // Gauge Link field GT is the gauge transform and lives on the VERTEX field
+  ////////////////////////////////////////////////////////////////////////////////////
+  void ForwardTriangles(GaugeField &Umu,LatticeComplex &plaq1,LatticeComplex &plaq2)
+  {
+    {
+      autoView(Umu_v,Umu,AcceleratorRead);
+      autoView(plaq1_v,plaq1,AcceleratorWrite);
+      autoView(plaq2_v,plaq2,AcceleratorWrite);
+      autoView(stencil_v,FaceStencil,AcceleratorRead);
+
+      accelerator_for(ss,EdgeGrid->oSites(),vComplex::Nsimd(),{
+
+	const int x = IcosahedronPatchX;
+	const int y = IcosahedronPatchY;
+	const int d = IcosahedronPatchDiagonal;
+	
+	auto Lx = Umu_v(ss)(x);
+	auto Ly = Umu_v(ss)(y);
+	auto Ld = Umu_v(ss)(d);
+
+	// for trace [ U_x(z) U_y(z+\hat x) adj(U_d(z)) ]
+	{
+	  auto SE1 = stencil_v.GetEntry(0,ss);
+	  auto doAdj = SE1->_adjoint;
+	  auto pol   = SE1->_polarisation;
+	  auto s1    = SE1->_offset;
+	  auto L1 = Umu_v(s1)(pol);
+	  if(doAdj)
+	    L1 = adj(L1);
+
+	  coalescedWrite(plaq1_v[ss](),trace(Lx*L1*adj(Ld) ) );
+	}
+	
+	// for trace [  U_y(z) adj(U_d(z))  U_x(z+\hat y) ]
+	{
+	  auto SE2 = stencil_v.GetEntry(1,ss);
+	  auto doAdj = SE2->_adjoint;
+	  auto pol   = SE2->_polarisation;
+	  auto s2    = SE2->_offset;
+	  auto L2 = Umu_v(s2)(pol);
+	  if(doAdj)
+	    L2 = adj(L2);
+
+	  coalescedWrite(plaq2_v[ss](),trace(Ly*adj(Ld)*L2 ) );
+	}
+      });
+    }
+  }
+
+  void  GaugeTransform(GaugeLinkField &gt, GaugeField &Umu)
+  {
+    assert(gt.Grid()==VertexGrid);
+    assert(Umu.Grid()==EdgeGrid);
+    
+    GridBase * vgrid = VertexGrid;
+    GridBase * grid  = EdgeGrid;
+
+    int osites  = grid->oSites();
+
+    uint64_t cart_sites = grid->CartesianOsites();
+    uint64_t Npole_sites = grid->NorthPoleOsites();
+    uint64_t Spole_sites = grid->SouthPoleOsites();
+    Coordinate pcoor = grid->ThisProcessorCoor();
+    Coordinate pgrid = grid->ProcessorGrid();
+    /*
+     * resize the stencil entries array and set npoints
+     */
+    autoView(g_v,gt,CpuRead);
+    autoView(Umu_v,Umu,CpuWrite);
+    for(uint64_t site=0;site<cart_sites; site ++) {
+
+      Coordinate Coor;
+      Coordinate NbrCoor;
+
+      int nd = grid->Nd();
+      int L  = grid->LocalDimensions()[0];
+	
+      ////////////////////////////////////////////////
+      // Outer index of neighbour Offset calculation
+      ////////////////////////////////////////////////
+      grid->oCoorFromOindex(Coor,site);
+      NbrCoor = Coor;
+      assert( grid->LocalDimensions()[1]==grid->LocalDimensions()[0]);
+      assert( grid->_simd_layout[0]==1); // Cannot vectorise in these dims
+      assert( grid->_simd_layout[1]==1);
+      assert( grid->_processors[0]==1); // Cannot mpi distribute in these dims
+      assert( grid->_processors[1]==1);
+
+      int Patch = Coor[nd-1];
+      int HemiPatch = Patch%HemiPatches;
+      int north = Patch/HemiPatches;
+      int south = 1-north;
+      int isPoleY;
+      int isPoleX;
+      
+      assert(Patch<IcosahedralPatches);
+      assert((north==1)||(south==1));
+
+      Coordinate XpCoor;
+      Coordinate YpCoor;
+      Coordinate DpCoor;
+
+      FaceStencil.GetNbrForPlusDiagonal(grid,Coor,DpCoor);
+      FaceStencil.GetNbrForPlusX(grid,Coor,XpCoor,isPoleX);
+      FaceStencil.GetNbrForPlusY(grid,Coor,YpCoor,isPoleY);
+
+      int XpHemiPatch  = XpCoor[nd-1]%HemiPatches;
+      int XpHemisphere = XpCoor[nd-1]/HemiPatches;
+
+      int DpPatch  = DpCoor[nd-1];
+      int DpHemiPatch  = DpCoor[nd-1]%HemiPatches;
+      int DpHemisphere = DpCoor[nd-1]/HemiPatches;
+
+      // Work out the pole_osite
+      Coordinate rdims;
+      Coordinate ocoor;
+      int64_t pole_osite;
+      int Ndm1 = grid->Nd()-1;
+      for(int d=2;d<Ndm1;d++){
+	int dd=d-2;
+	rdims.push_back(grid->_rdimensions[d]);
+	ocoor.push_back(Coor[d]%grid->_rdimensions[d]);
+      }
+      Lexicographic::IndexFromCoor(ocoor,pole_osite,rdims);
+      
+      uint64_t xp_idx;
+      uint64_t yp_idx;
+      uint64_t dp_idx;
+      if ( isPoleX ) {
+	assert(vgrid->ownsSouthPole());
+	xp_idx = pole_osite + grid->SouthPoleOsite();
+      } else {
+	xp_idx = grid->oIndex(XpCoor);
+      }
+      if ( isPoleY ) {
+	assert(vgrid->ownsNorthPole());
+	yp_idx = pole_osite + grid->NorthPoleOsite();
+      } else {
+	yp_idx = grid->oIndex(YpCoor);
+      }
+      dp_idx = grid->oIndex(DpCoor);
+
+      auto g  = g_v(site)();
+      auto gx = g_v(xp_idx)();
+      auto gy = g_v(yp_idx)();
+      auto gd = g_v(dp_idx)();
+
+      auto lx = Umu_v(site)(IcosahedronPatchX);
+      auto ly = Umu_v(site)(IcosahedronPatchY);
+      auto ld = Umu_v(site)(IcosahedronPatchDiagonal);
+      
+      lx = g*lx*adj(gx);
+      ly = g*ly*adj(gy);
+      ld = g*ld*adj(gd);
+
+      coalescedWrite(Umu_v[site](IcosahedronPatchX),lx);
+      coalescedWrite(Umu_v[site](IcosahedronPatchY),ly);
+      coalescedWrite(Umu_v[site](IcosahedronPatchDiagonal),ld);
+    };
+  }
+
+};
 
 int main (int argc, char ** argv)
 {
   Grid_init(&argc,&argv);
 
   const int L=8;
-  const int Npatch = num_icosahedron_tiles;
+  const int Npatch = IcosahedralPatches;
 
   // Put SIMD all in time direction
   Coordinate latt_size = GridDefaultLatt();
-  Coordinate simd_layout({1,1,vComplexD::Nsimd(),1});
+  Coordinate simd_layout({1,1,vComplex::Nsimd(),1});
   Coordinate mpi_layout = GridDefaultMpi();
 
   std::cout << GridLogMessage << " mpi "<<mpi_layout<<std::endl;
@@ -57,7 +261,7 @@ int main (int argc, char ** argv)
   GridCartesianCrossIcosahedron VertexGrid(latt_size,simd_layout,mpi_layout,IcosahedralVertices);
 
   std::cout << GridLogMessage << " Created vertex grid "<<std::endl;
-  LatticeIcosahedralLorentzComplexD Umu(&EdgeGrid);
+  LatticeIcosahedralLorentzColourMatrix Umu(&EdgeGrid);
   LatticeComplex Phi(&VertexGrid);
   std::cout << GridLogMessage << " Created two fields "<<std::endl;
 
@@ -65,7 +269,7 @@ int main (int argc, char ** argv)
   Umu = Zero();
   std::cout << GridLogMessage << " Zeroed two fields "<<std::endl;
 
-  ComplexD one (1.0);
+  Complex one (1.0);
   Phi = one;
   Umu = one;
   
@@ -78,16 +282,55 @@ int main (int argc, char ** argv)
   //  std::cout << " Umu "<<Umu<<std::endl;
   //  std::cout << " Phi "<<Phi<<std::endl;
   LatticePole(Phi,South);
-  std::cout << " Phi South Pole set\n"<<Phi<<std::endl;
+  //  std::cout << " Phi South Pole set\n"<<Phi<<std::endl;
 
   LatticePole(Phi,North);
-  std::cout << " Phi North Pole set\n"<<Phi<<std::endl;
+  //  std::cout << " Phi North Pole set\n"<<Phi<<std::endl;
 
   for(int mu=0;mu<VertexGrid._ndimension;mu++){
     std::cout << " Calling lattice coordinate mu="<<mu<<std::endl;
     LatticeCoordinate(Phi,mu);
-    std::cout << " Phi coor mu="<<mu<<"\n"<<Phi<<std::endl;
+    //    std::cout << " Phi coor mu="<<mu<<"\n"<<Phi<<std::endl;
   }
 
+  std::cout << "Creating face stencil"<<std::endl;
+  IcosahedralEdgeSupport<IcosahedralGimpl> Support(&VertexGrid,&EdgeGrid);
+
+  std::cout << " Calling Test Geometry "<<std::endl;
+  Support.FaceStencil.TestGeometry();
+
+
+  Umu=one;
+  LatticeComplex plaq1(&EdgeGrid);
+  LatticeComplex plaq2(&EdgeGrid);
+
+  Support.ForwardTriangles(Umu,plaq1,plaq2);
+  std::cout << " plaq1 "<< norm2(plaq1)<<std::endl;
+  std::cout << " plaq2 "<< norm2(plaq2)<<std::endl;
+
+  // Random gauge xform
+  std::vector<int> seeds({1,2,3,4});
+  GridParallelRNG  vRNG(&EdgeGrid);   vRNG.SeedFixedIntegers(seeds);
+  LatticeIcosahedralColourMatrix g(&VertexGrid);
+  
+  //  SU<Nc>::LieRandomize(vRNG,g);
+
+  LatticeReal gr(&VertexGrid);
+  LatticeComplex gc(&VertexGrid);
+  gaussian(vRNG,gr);
+  Complex ci(0.0,1.0);
+  gc = toComplex(gr);
+  g=one;
+  g = g * exp(ci*gc);
+  
+  std::cout << "applying gauge transform"<<std::endl;
+  Support.GaugeTransform(g,Umu);
+  std::cout << "applied gauge transform "<<Umu<<std::endl;
+  
+  Support.ForwardTriangles(Umu,plaq1,plaq2);
+  std::cout << " plaq1 "<< norm2(plaq1)<<std::endl;
+  std::cout << " plaq2 "<< norm2(plaq2)<<std::endl;
+
+  
   Grid_finalize();
 }