From f9cda24781e6f7b3291947e7ae6b4597daa91375 Mon Sep 17 00:00:00 2001
From: Quadro <Peter Boyle>
Date: Tue, 8 Jun 2021 21:50:43 -0400
Subject: [PATCH] Cleaner

---
 Grid/qcd/action/momentum/DDHMCfilter.h | 62 ++------------------------
 1 file changed, 4 insertions(+), 58 deletions(-)

diff --git a/Grid/qcd/action/momentum/DDHMCfilter.h b/Grid/qcd/action/momentum/DDHMCfilter.h
index 716cce1d..973a3267 100644
--- a/Grid/qcd/action/momentum/DDHMCfilter.h
+++ b/Grid/qcd/action/momentum/DDHMCfilter.h
@@ -36,71 +36,17 @@ NAMESPACE_BEGIN(Grid);
 template<typename MomentaField>
 struct DDHMCFilter: public MomentumFilterBase<MomentaField>
 {
-  //  typedef typename MomentaField::vector_type vector_type; //SIMD-vectorized complex type
-  //  typedef typename MomentaField::scalar_type scalar_type; //scalar complex type
-  //  typedef iVector<iScalar<iScalar<vector_type> >, Nd > LorentzScalarType; //complex phase for each site/direction
-  //  typedef iScalar<iScalar<iScalar<vector_type> > >            ScalarType; //complex phase for each site
-  //  typedef Lattice<LorentzScalarType> LatticeLorentzScalarType;
-  //  typedef Lattice<ScalarType> LatticeScalarType;
-  
   Coordinate Block;
   int Width;
-
-  // Could also zero links in plane like Luscher advocates.
   
-  DDHMCFilter(const Coordinate &_Block,int _Width): Block(_Block), Width(_Width)
-  {
-    assert( (Width==0) || (Width==1) || (Width==2) );
-  }
+  DDHMCFilter(const Coordinate &_Block): Block(_Block) {}
 
   void applyFilter(MomentaField &P) const override
   {
-    GridBase *grid = P.Grid();
-
-    LatticeColourMatrix zz(grid); zz = Zero();
-    MomentaField zzz(grid); zzz = Zero();
-
-    ////////////////////////////////////////////////////
-    // Zero strictly links crossing between domains
-    // Luscher also zeroes links in plane of domain boundaries
-    // Keeping interior only. This prevents force from plaquettes
-    // crossing domains and keeps whole MD trajectory local.
-    // This is the case Width=1.
-    // Width = 2 should also decouple rectangle force.
-    ////////////////////////////////////////////////////
-    Coordinate Global=grid->GlobalDimensions();
-    LatticeInteger coor(grid); 
-    for(int mu=0;mu<Nd;mu++) {
-
-#define PROJECT_DOMAIN      
-      if ( (Block[mu] <= Global[mu])&&(Block[mu]>0) ) {
-	LatticeCoordinate(coor,mu);
-#if 1
-	P = where(mod(coor,Block[mu])==Integer(Block[mu]-3),zzz,P); //width 4
-	P = where(mod(coor,Block[mu])==Integer(Block[mu]-2),zzz,P); //width 4
-	P = where(mod(coor,Block[mu])==Integer(Block[mu]-1),zzz,P); //width 2
-	P = where(mod(coor,Block[mu])==Integer(0),zzz,P);           //width 2
-	P = where(mod(coor,Block[mu])==Integer(1),zzz,P);           //width 4
-	auto P_mu = PeekIndex<LorentzIndex>(P,mu);                  
-	P = where(mod(coor,Block[mu])==Integer(2),zzz,P);           //width 6
-	PokeIndex<LorentzIndex>(P, P_mu, mu);
-#else
-	P = where(mod(coor,Block[mu])==Integer(Block[mu]-1),zzz,P); //width 2
-	auto P_mu = PeekIndex<LorentzIndex>(P,mu);                  
-	P = where(mod(coor,Block[mu])==Integer(0),zzz,P);           //width 6
-	PokeIndex<LorentzIndex>(P, P_mu, mu);
-#endif	
-      }      
-    }
-#ifdef PROJECT_DOMAIN    
-    LatticeInteger domaincb(grid); domaincb=Zero();
-    for(int d=0;d<Nd;d++){
-      LatticeCoordinate(coor,d);
-      domaincb = domaincb + div(coor,Block[d]);
-    }
-    P = where(mod(domaincb,2)==Integer(1),P,zzz);
-#endif    
+    DomainDecomposition Domains(Block);
+    Domains.ProjectDDHMC(P);
   }
 };
+
 NAMESPACE_END(Grid);