Fixing Bug in DCPSE

src/DCPSE_op/DCPSE_Solver.hpp

@@ -495,7 +495,7 @@ public:
     typename Sys_eqs::Vector_type & getB(options_solver opt = options_solver::STANDARD)
     {
 #ifdef SE_CLASS1
-        consistency();
+        //consistency();
 #endif
         if (opt == options_solver::LAGRANGE_MULTIPLIER)
         {

src/DCPSE_op/DCPSE_copy/DcpseRhs.hpp

@@ -51,10 +51,10 @@ MatrixType &DcpseRhs<dim>::getVector(MatrixType &b)
         b(i, 0) = sign * dm.evaluate(Point<dim, T>(0));
     }
     //Choosing a(0,0) as  a free parameter can let us set b(0,0) for numerical robustness
-    if (b(0,0) == 0.0 && sign == 1)
-    {
-        b(0,0) = 100;
-    }
+    //if (b(0,0) == 0.0 && sign == 1)
+    //{
+    //    b(0,0) = 100;
+    //}
 
 
 

src/DCPSE_op/DCPSE_op.hpp

@@ -405,10 +405,10 @@ public:
                 o1.template value_nz<Sys_eqs>(p_map,k,cols,coeff_k,comp);
 
                 auto coeff_kk = dcp[i].getSign() * coeff_dc * coeff / dcp[i].getEpsilonPrefactor(key);
-                o1.template value_nz<Sys_eqs>(p_map,key,cols,coeff_k,comp);
+                o1.template value_nz<Sys_eqs>(p_map,key,cols,coeff_kk,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>(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);*/
             }
         }
     }
@@ -583,7 +583,7 @@ public:
                 o2.template value_nz<Sys_eqs>(p_map,k,cols,k_coeff,comp);
 
                 auto kk_coeff = o1.value(key) * dcp[i].getSign() * coeff_dc * coeff / dcp[i].getEpsilonPrefactor(key);
-                o2.template value_nz<Sys_eqs>(p_map,k,cols,k_coeff,comp);
+                o2.template value_nz<Sys_eqs>(p_map,k,cols,kk_coeff,comp);
 
                 //cols[p_map. template getProp<0>(k)*Sys_eqs::nvar + comp] += o1.value(key)[i] * coeff_dc * coeff / dcp[i].getEpsilonPrefactor(key);
                 //cols[p_map. template getProp<0>(key)*Sys_eqs::nvar + comp] += o1.value(key)[i] * dcp[i].getSign() * coeff_dc * coeff / dcp[i].getEpsilonPrefactor(key);