Fixing conflicts

src/DCPSE_op/DCPSE_Solver.hpp

@@ -224,12 +224,12 @@ class DCPSE_scheme: public MatrixAssembler
         	std::cerr << "Error " << __FILE__ << ":" << __LINE__ << " the term B and the Matrix A for Ax=B must contain the same number of rows\n";
         	return;
         }
-
-        if (row_b != p_map.size_local() * Sys_eqs::nvar)
-        {
-        	std::cerr << "Error " << __FILE__ << ":" << __LINE__ << " your system is underdetermined you set " << row_b << " conditions " << " but i am expecting " << p_map.size_local() * Sys_eqs::nvar << std::endl;
-        	return;
-        }
+            if (row_b != p_map.size_local() * Sys_eqs::nvar) {
+                std::cerr << "Error " << __FILE__ << ":" << __LINE__ << " your system is underdetermined you set "
+                          << row_b << " conditions " << " but i am expecting " << p_map.size_local() * Sys_eqs::nvar
+                          << std::endl;
+                return;
+            }
 
         // Indicate all the non zero rows
         openfpm::vector<unsigned char> nz_rows;
@@ -443,7 +443,6 @@ public:
             A.resize(tot * Sys_eqs::nvar, tot * Sys_eqs::nvar,
                      p_map.size_local() * Sys_eqs::nvar,
                      p_map.size_local() * Sys_eqs::nvar);
-//        std::cout << Eigen::MatrixXd(A) << std::endl;
         } else
         {
             auto & v_cl = create_vcluster();
@@ -483,6 +482,7 @@ public:
 
 
         }
+
         return A;
 
     }
@@ -506,6 +506,7 @@ public:
             }
         }
 
+
         return b;
     }
 

src/DCPSE_op/DCPSE_copy/MonomialBasis.hpp

@@ -124,9 +124,9 @@ void MonomialBasis<dim>::generateBasis(std::vector<unsigned int> m, unsigned int
     grid_key_dx_iterator_sub_bc<dim> it(grid, start, stop, bc);
 
     // Finally compute alpha_min
-    unsigned char alphaMin = static_cast<unsigned char>(!(mSum % 2)); // if mSum is even, alpha_min must be 1
+    //unsigned char alphaMin = static_cast<unsigned char>(!(mSum % 2)); // if mSum is even, alpha_min must be 1
     //std::cout<<"AlphaMin: "<<alphaMin<<std::endl;
-    //unsigned char alphaMin = 0; // we want to always have 1 in the basis
+    unsigned char alphaMin = 0; // we want to always have 1 in the basis
 
     while (it.isNext())
     {

src/DCPSE_op/DCPSE_op.hpp

@@ -65,9 +65,15 @@ public:
             {
                 auto coeff_dc = dcp.getCoeffNN(key,j);
                 auto k = dcp.getIndexNN(key,j);
-                cols[p_map. template getProp<0>(k)*Sys_eqs::nvar + comp] += coeff_dc * coeff / dcp.getEpsilonPrefactor(key);
 
-                cols[p_map. template getProp<0>(key)*Sys_eqs::nvar + comp] += dcp.getSign() * coeff_dc * coeff / dcp.getEpsilonPrefactor(key);
+                auto k_coeff = coeff_dc * coeff / dcp[i].getEpsilonPrefactor(key);
+                o1.template value_nz<Sys_eqs>(p_map,k,cols,k_coeff,comp);
+
+                auto kk_coeff = dcp[i].getSign() * coeff_dc * coeff / dcp[i].getEpsilonPrefactor(key);
+                o1.template value_nz<Sys_eqs>(p_map,k,cols,k_coeff,comp);
+                //cols[p_map. template getProp<0>(k)*Sys_eqs::nvar + comp] += coeff_dc * coeff / dcp.getEpsilonPrefactor(key);
+
+                //cols[p_map. template getProp<0>(key)*Sys_eqs::nvar + comp] += dcp.getSign() * coeff_dc * coeff / dcp.getEpsilonPrefactor(key);
             }
     }
 
@@ -156,10 +162,16 @@ public:
                 auto coeff_dc = dcp[i].getCoeffNN(key,j);
                 auto k = dcp[i].getIndexNN(key,j);
 
+                auto k_coeff = coeff_dc * coeff / dcp[i].getEpsilonPrefactor(key);
+                o1.template value_nz<Sys_eqs>(p_map,k,cols,k_coeff,comp);
+
+                auto kk_coeff = dcp[i].getSign() * coeff_dc * coeff / dcp[i].getEpsilonPrefactor(key);
+                o1.template value_nz<Sys_eqs>(p_map,k,cols,k_coeff,comp);
+
 
-                cols[p_map. template getProp<0>(k)*Sys_eqs::nvar + comp] += coeff_dc * coeff / dcp[i].getEpsilonPrefactor(key);
+                //cols[p_map. template getProp<0>(k)*Sys_eqs::nvar + comp] += coeff_dc * coeff / dcp[i].getEpsilonPrefactor(key);
 
-                cols[p_map. template getProp<0>(key)*Sys_eqs::nvar + comp] += dcp[i].getSign() * coeff_dc * coeff / dcp[i].getEpsilonPrefactor(key);
+                //cols[p_map. template getProp<0>(key)*Sys_eqs::nvar + comp] += dcp[i].getSign() * coeff_dc * coeff / dcp[i].getEpsilonPrefactor(key);
             }
         }
     }
@@ -249,10 +261,16 @@ public:
                 auto coeff_dc = dcp[i].getCoeffNN(key,j);
                 auto k = dcp[i].getIndexNN(key,j);
 
+                auto k_coeff = coeff_dc * coeff / dcp[i].getEpsilonPrefactor(key);
+                o1.template value_nz<Sys_eqs>(p_map,k,cols,k_coeff,comp);
 
-                cols[p_map. template getProp<0>(k)*Sys_eqs::nvar + comp] += coeff_dc * coeff / dcp[i].getEpsilonPrefactor(key);
+                auto kk_coeff = dcp[i].getSign() * coeff_dc * coeff / dcp[i].getEpsilonPrefactor(key);
+                o1.template value_nz<Sys_eqs>(p_map,k,cols,k_coeff,comp);
 
-                cols[p_map. template getProp<0>(key)*Sys_eqs::nvar + comp] += dcp[i].getSign() * coeff_dc * coeff / dcp[i].getEpsilonPrefactor(key);
+
+                //cols[p_map. template getProp<0>(k)*Sys_eqs::nvar + comp] += coeff_dc * coeff / dcp[i].getEpsilonPrefactor(key);
+
+                //cols[p_map. template getProp<0>(key)*Sys_eqs::nvar + comp] += dcp[i].getSign() * coeff_dc * coeff / dcp[i].getEpsilonPrefactor(key);
             }
         }
     }
@@ -490,10 +508,16 @@ public:
                 auto coeff_dc = dcp[i].getCoeffNN(key,j);
                 auto k = dcp[i].getIndexNN(key,j);
 
+                auto k_coeff = coeff_dc * coeff / dcp[i].getEpsilonPrefactor(key);
+                o1.template value_nz<Sys_eqs>(p_map,k,cols,k_coeff,comp);
+
+                auto kk_coeff = dcp[i].getSign() * coeff_dc * coeff / dcp[i].getEpsilonPrefactor(key);
+                o1.template value_nz<Sys_eqs>(p_map,k,cols,k_coeff,comp);
+
 
-                cols[p_map. template getProp<0>(k)*Sys_eqs::nvar + comp] += coeff_dc * coeff / dcp[i].getEpsilonPrefactor(key);
+                //cols[p_map. template getProp<0>(k)*Sys_eqs::nvar + comp] += coeff_dc * coeff / dcp[i].getEpsilonPrefactor(key);
 
-                cols[p_map. template getProp<0>(key)*Sys_eqs::nvar + comp] += dcp[i].getSign() * coeff_dc * coeff / dcp[i].getEpsilonPrefactor(key);
+               //cols[p_map. template getProp<0>(key)*Sys_eqs::nvar + comp] += dcp[i].getSign() * coeff_dc * coeff / dcp[i].getEpsilonPrefactor(key);
             }
         }
     }

src/DCPSE_op/DCPSE_op_test2.cpp

@@ -188,15 +188,15 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests2)
         Divergence Div(Particles, 2, rCut, 1.9);
         int n=10;
         double nu=1e-2;
-        for(int i=0; i<=n ;i++)
+        Particles.write_frame("Stokes",0);
+        for(int i=1; i<=n ;i++)
         {   RHS=-Grad(P);
             DCPSE_scheme<equations2,decltype(Particles)> Solver(2 * rCut, Particles);
-//            auto Stokes = Adv(V,V_star)-nu*Lap(V_star);
             auto Stokes1 = Adv(V[0],V_star[0])-nu*Lap(V_star[0]);
             auto Stokes2 = Adv(V[1],V_star[1])-nu*Lap(V_star[1]);
             Solver.impose(Stokes1,bulk,RHS[0],vx);
             Solver.impose(Stokes2,bulk,RHS[1],vy);
-            Solver.impose(V_star[0], up_p,1,vx);
+            Solver.impose(V_star[0], up_p,1.0,vx);
             Solver.impose(V_star[1], up_p,0,vy);
             Solver.impose(V_star[0], r_p, 0,vx);
             Solver.impose(V_star[1], r_p, 0,vy);