From 5617e3d5a35c8d4de08dc3637f971d93e9a7205f Mon Sep 17 00:00:00 2001
From: Pietro Incardona <incardon@mpi-cbg.de>
Date: Thu, 26 Mar 2020 17:31:57 +0100
Subject: [PATCH] DCPSE solver working with system of equations

---
 src/DCPSE_op/DCPSE_Solver.hpp                 | 57 +++++++++++--
 src/DCPSE_op/DCPSE_op.hpp                     | 64 +++++++--------
 src/DCPSE_op/DCPSE_op_test2.cpp               | 25 ++++--
 .../Vector/vector_dist_operators.hpp          | 82 ++++++++++++-------
 4 files changed, 156 insertions(+), 72 deletions(-)

diff --git a/src/DCPSE_op/DCPSE_Solver.hpp b/src/DCPSE_op/DCPSE_Solver.hpp
index 64be0493..f6cc81de 100644
--- a/src/DCPSE_op/DCPSE_Solver.hpp
+++ b/src/DCPSE_op/DCPSE_Solver.hpp
@@ -16,6 +16,27 @@
 template<unsigned int prp_id>
 struct prop_id {};
 
+class eq_id
+{
+	int id;
+
+public:
+
+	eq_id()
+	:id(0)
+	{}
+
+	int getId()
+	{
+		return id;
+	}
+
+	void setId(int id)
+	{
+		this->id = id;
+	}
+};
+
 enum options_solver
 {
     STANDARD,
@@ -189,6 +210,7 @@ class DCPSE_scheme: public MatrixAssembler
         }
     };
 
+
     /*! \brief Check if the Matrix is consistent
  *
  */
@@ -327,13 +349,38 @@ public:
     template<typename T, typename index_type, unsigned int prp_id>
     void impose(const T & op , openfpm::vector<index_type> & subset,
                                                           const prop_id<prp_id> & num,
-                                                          long int id = 0)
+                                                          eq_id id = eq_id())
     {
         auto itd = subset.template getIteratorElements<0>();
 
         variable_b<prp_id> vb(parts);
 
-        impose_git(op,vb,id,itd);
+        impose_git(op,vb,id.getId(),itd);
+    }
+
+    /*! \brief Impose an operator
+*
+* This function impose an operator on a particular grid region to produce the system
+*
+* Ax = b
+*
+* ## Stokes equation 2D, lid driven cavity with one splipping wall
+* \snippet eq_unit_test.hpp Copy the solution to grid
+*
+* \param op Operator to impose (A term)
+* \param num right hand side of the term (b term)
+* \param id Equation id in the system that we are imposing
+* \param it_d iterator that define where you want to impose
+*
+*/
+    template<typename T, typename index_type, typename RHS_type, typename sfinae = typename std::enable_if< !std::is_fundamental<RHS_type>::type::value >::type >
+    void impose(const T & op , openfpm::vector<index_type> & subset,
+                                                          const RHS_type & rhs,
+                                                          eq_id id = eq_id())
+    {
+        auto itd = subset.template getIteratorElements<0>();
+
+        impose_git(op,rhs,id.getId(),itd);
     }
 
     /*! \brief Impose an operator
@@ -354,13 +401,13 @@ public:
     template<typename T, typename index_type> void impose(const T & op ,
                                                                           openfpm::vector<index_type> & subset,
                                                                           const typename Sys_eqs::stype num,
-                                                                          long int id = 0)
+                                                                          eq_id id = eq_id())
     {
         auto itd = subset.template getIteratorElements<0>();
 
         constant_b b(num);
 
-        impose_git(op,b,id,itd);
+        impose_git(op,b,id.getId(),itd);
     }
 
     /*! \brief produce the Matrix
@@ -477,7 +524,7 @@ public:
 
             // Calculate the non-zero colums
             typename Sys_eqs::stype coeff = 1.0;
-            op.value_nz(p_map,key,cols,coeff);
+            op.template value_nz<Sys_eqs>(p_map,key,cols,coeff,0);
 
             // indicate if the diagonal has been set
             bool is_diag = false;
diff --git a/src/DCPSE_op/DCPSE_op.hpp b/src/DCPSE_op/DCPSE_op.hpp
index 8e2d7699..abe53792 100644
--- a/src/DCPSE_op/DCPSE_op.hpp
+++ b/src/DCPSE_op/DCPSE_op.hpp
@@ -57,17 +57,17 @@ public:
         return dcp.computeDifferentialOperator(key,o1);
     }
 
-    template<typename pmap_type, typename unordered_map_type, typename coeff_type>
-    inline void value_nz(pmap_type & p_map, const vect_dist_key_dx & key, unordered_map_type & cols, coeff_type & coeff) const
+    template<typename Sys_eqs, typename pmap_type, typename unordered_map_type, typename coeff_type>
+    inline void value_nz(pmap_type & p_map, const vect_dist_key_dx & key, unordered_map_type & cols, coeff_type & coeff, unsigned int comp) const
     {
             // for all NN of key
             for (int j = 0 ; j < dcp.getNumNN(key) ; j++)
             {
                 auto coeff_dc = dcp.getCoeffNN(key,j);
                 auto k = dcp.getIndexNN(key,j);
-                cols[p_map. template getProp<0>(k)] += coeff_dc * coeff / dcp.getEpsilonPrefactor(key);
+                cols[p_map. template getProp<0>(k)*Sys_eqs::nvar + comp] += coeff_dc * coeff / dcp.getEpsilonPrefactor(key);
 
-                cols[p_map. template getProp<0>(key)] += dcp.getSign() * 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);
             }
     }
 
@@ -145,8 +145,8 @@ public:
         return v_grad;
     }
 
-    template<typename pmap_type, typename unordered_map_type, typename coeff_type>
-    inline void value_nz(pmap_type & p_map, const vect_dist_key_dx & key, unordered_map_type & cols, coeff_type & coeff) const
+    template<typename Sys_eqs, typename pmap_type, typename unordered_map_type, typename coeff_type>
+    inline void value_nz(pmap_type & p_map, const vect_dist_key_dx & key, unordered_map_type & cols, coeff_type & coeff, unsigned int comp) const
     {
         for (int i = 0 ; i < DCPSE_type::vtype::dims ; i++)
         {
@@ -157,9 +157,9 @@ public:
                 auto k = dcp[i].getIndexNN(key,j);
 
 
-                cols[p_map. template getProp<0>(k)] += 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)] += 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);
             }
         }
     }
@@ -238,8 +238,8 @@ public:
         return v_grad;
     }
 
-    template<typename pmap_type, typename unordered_map_type, typename coeff_type>
-    inline void value_nz(pmap_type & p_map, const vect_dist_key_dx & key, unordered_map_type & cols, coeff_type & coeff) const
+    template<typename Sys_eqs, typename pmap_type, typename unordered_map_type, typename coeff_type>
+    inline void value_nz(pmap_type & p_map, const vect_dist_key_dx & key, unordered_map_type & cols, coeff_type & coeff, unsigned int comp) const
     {
         for (int i = 0 ; i < DCPSE_type::vtype::dims ; i++)
         {
@@ -250,9 +250,9 @@ public:
                 auto k = dcp[i].getIndexNN(key,j);
 
 
-                cols[p_map. template getProp<0>(k)] += 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)] += 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);
             }
         }
     }
@@ -300,17 +300,17 @@ public:
             :l1(l1),l2(l2)
     {}
 
-    template<typename pmap_type, typename unordered_map_type, typename coeff_type>
-    inline void value_nz(pmap_type & p_map, const vect_dist_key_dx & key, unordered_map_type & cols, coeff_type & coeff) const
+    template<typename Sys_eqs, typename pmap_type, typename unordered_map_type, typename coeff_type>
+    inline void value_nz(pmap_type & p_map, const vect_dist_key_dx & key, unordered_map_type & cols, coeff_type & coeff, unsigned int comp) const
     {
         if (l1.template get<0>(i) != key.getKey())
         {
             std::cout << "ERROR" << std::endl;
         }
 
-        cols[p_map. template getProp<0>(key)] += coeff;
+        cols[p_map. template getProp<0>(key)*Sys_eqs::nvar + comp] += coeff;
         std::cout << "L2: " << l2.template get<0>(i) << std::endl;
-        cols[p_map. template getProp<0>(l2.template get<0>(i))] -= coeff;
+        cols[p_map. template getProp<0>(l2.template get<0>(i))*Sys_eqs::nvar + comp ] -= coeff;
 
         i++;
     }
@@ -333,16 +333,16 @@ public:
     :l1(l1),l2_key(l2_key)
     {}
 
-    template<typename pmap_type, typename unordered_map_type, typename coeff_type>
-    inline void value_nz(pmap_type & p_map, const vect_dist_key_dx & key, unordered_map_type & cols, coeff_type & coeff)  const
+    template<typename Sys_eqs, typename pmap_type, typename unordered_map_type, typename coeff_type>
+    inline void value_nz(pmap_type & p_map, const vect_dist_key_dx & key, unordered_map_type & cols, coeff_type & coeff, unsigned int comp)  const
     {
         if (l1.template get<0>(i) != key.getKey())
         {
             std::cout << "ERROR" << std::endl;
         }
 
-        cols[p_map. template getProp<0>(key)] += coeff;
-        cols[p_map. template getProp<0>(l2_key)] -= coeff;
+        cols[p_map. template getProp<0>(key)*Sys_eqs::nvar + comp] += coeff;
+        cols[p_map. template getProp<0>(l2_key)*Sys_eqs::nvar + comp] -= coeff;
         i++;
     }
 };
@@ -390,8 +390,8 @@ public:
         return v_lap;
     }
 
-    template<typename pmap_type, typename unordered_map_type, typename coeff_type>
-    inline void value_nz(pmap_type & p_map, const vect_dist_key_dx & key, unordered_map_type & cols, coeff_type & coeff) const
+    template<typename Sys_eqs, typename pmap_type, typename unordered_map_type, typename coeff_type>
+    inline void value_nz(pmap_type & p_map, const vect_dist_key_dx & key, unordered_map_type & cols, coeff_type & coeff, unsigned int comp) const
     {
         for (int i = 0 ; i < DCPSE_type::vtype::dims ; i++)
         {
@@ -402,9 +402,9 @@ public:
                 auto k = dcp[i].getIndexNN(key,j);
 
 
-                cols[p_map. template getProp<0>(k)] += 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)] += 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);
             }
         }
     }
@@ -475,8 +475,8 @@ public:
         return v_div;
     }
 
-    template<typename pmap_type, typename unordered_map_type, typename coeff_type>
-    inline void value_nz(pmap_type & p_map, const vect_dist_key_dx & key, unordered_map_type & cols, coeff_type & coeff) const
+    template<typename Sys_eqs, typename pmap_type, typename unordered_map_type, typename coeff_type>
+    inline void value_nz(pmap_type & p_map, const vect_dist_key_dx & key, unordered_map_type & cols, coeff_type & coeff, unsigned int comp) const
     {
         for (int i = 0 ; i < DCPSE_type::vtype::dims ; i++)
         {
@@ -487,9 +487,9 @@ public:
                 auto k = dcp[i].getIndexNN(key,j);
 
 
-                cols[p_map. template getProp<0>(k)] += 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)] += 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);
             }
         }
     }
@@ -564,8 +564,8 @@ public:
     }
 
 
-    template<typename pmap_type, typename unordered_map_type, typename coeff_type>
-    inline void value_nz(pmap_type & p_map, const vect_dist_key_dx & key, unordered_map_type & cols, coeff_type & coeff) const
+    template<typename Sys_eqs, typename pmap_type, typename unordered_map_type, typename coeff_type>
+    inline void value_nz(pmap_type & p_map, const vect_dist_key_dx & key, unordered_map_type & cols, coeff_type & coeff, unsigned int comp) const
     {
         for (int i = 0 ; i < DCPSE_type::vtype::dims ; i++)
         {
@@ -576,9 +576,9 @@ public:
                 auto k = dcp[i].getIndexNN(key,j);
 
 
-                cols[p_map. template getProp<0>(k)] += 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)] += 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);
             }
         }
     }
diff --git a/src/DCPSE_op/DCPSE_op_test2.cpp b/src/DCPSE_op/DCPSE_op_test2.cpp
index 761a9229..34508295 100644
--- a/src/DCPSE_op/DCPSE_op_test2.cpp
+++ b/src/DCPSE_op/DCPSE_op_test2.cpp
@@ -26,7 +26,7 @@ struct equations2 {
     //! dimensionaly of the equation ( 3D problem ...)
     static const unsigned int dims = 2;
     //! number of fields in the system
-    static const unsigned int nvar = 1;
+    static const unsigned int nvar = 2;
 
     //! boundary at X and Y
     static const bool boundary[];
@@ -150,6 +150,12 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests2)
         }
         V_BC=V;
 
+
+        eq_id vx,vy;
+
+        vx.setId(0);
+        vy.setId(1);
+
         Derivative_x Dx(Particles, 2, rCut,1.9);
         Derivative_y Dy(Particles, 2, rCut,1.9);
         Gradient Grad(Particles, 2, rCut,1.9);
@@ -161,12 +167,17 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests2)
         for(int i=0; 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);
-            Solver.impose(Stokes,bulk,prop_id<3>());
-            Solver.impose(V_star, up_p,prop_id<4>());
-            Solver.impose(V_star, r_p, 0);
-            Solver.impose(V_star, dw_p,0);
-            Solver.impose(V_star, l_p, 0);
+//            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[1], up_p,0,vy);
+            Solver.impose(V_star[0], r_p, 0,vx);
+            Solver.impose(V_star[1], r_p, 0,vy);
+            Solver.impose(V_star, dw_p,0,vx);
+            Solver.impose(V_star, l_p, 0,vy);
             Solver.solve(V_star);
             std::cout << "Stokes Solved" << std::endl;
             RHS=Div(V_star);
diff --git a/src/Operators/Vector/vector_dist_operators.hpp b/src/Operators/Vector/vector_dist_operators.hpp
index 09798675..64ee1efc 100644
--- a/src/Operators/Vector/vector_dist_operators.hpp
+++ b/src/Operators/Vector/vector_dist_operators.hpp
@@ -231,11 +231,11 @@ public:
 		return o1.value(key) + o2.value(key);
 	}
 
-    template<typename pmap_type, typename unordered_map_type, typename coeff_type>
-    inline void value_nz(pmap_type & p_map, const vect_dist_key_dx & key, unordered_map_type & cols, coeff_type & coeff) const
+    template<typename Sys_eqs, typename pmap_type, typename unordered_map_type, typename coeff_type>
+    inline void value_nz(pmap_type & p_map, const vect_dist_key_dx & key, unordered_map_type & cols, coeff_type & coeff, unsigned int comp) const
     {
-        o1.value_nz(p_map,key,cols,coeff);
-        o2.value_nz(p_map,key,cols,coeff);
+        o1.template value_nz<Sys_eqs>(p_map,key,cols,coeff, comp);
+        o2.template value_nz<Sys_eqs>(p_map,key,cols,coeff, comp);
     }
     /*! \brief Return the vector on which is acting
      *
@@ -334,12 +334,12 @@ public:
         return first_or_second<has_vtype<exp1>::value,exp1,exp2>::getVector(o1,o2);
     }
 
-    template<typename pmap_type, typename unordered_map_type, typename coeff_type>
-    inline void value_nz(pmap_type & p_map, const vect_dist_key_dx & key, unordered_map_type & cols, coeff_type & coeff) const
+    template<typename Sys_eqs,typename pmap_type, typename unordered_map_type, typename coeff_type>
+    inline void value_nz(pmap_type & p_map, const vect_dist_key_dx & key, unordered_map_type & cols, coeff_type & coeff, unsigned int comp) const
     {
-        o1.value_nz(p_map,key,cols,coeff);
+        o1.template value_nz<Sys_eqs>(p_map,key,cols,coeff,comp);
         coeff_type tmp = -coeff;
-        o2.value_nz(p_map,key,cols,tmp);
+        o2.template value_nz<Sys_eqs>(p_map,key,cols,tmp,comp);
     }
 };
 
@@ -391,11 +391,11 @@ public:
 		return o1.value(key) * o2.value(key);
 	}
 
-    template<typename pmap_type, typename unordered_map_type, typename coeff_type>
-    inline void value_nz(pmap_type & p_map, const vect_dist_key_dx & key, unordered_map_type & cols, coeff_type & coeff) const
+    template<typename Sys_eqs,typename pmap_type, typename unordered_map_type, typename coeff_type>
+    inline void value_nz(pmap_type & p_map, const vect_dist_key_dx & key, unordered_map_type & cols, coeff_type & coeff, unsigned int comp) const
     {
-        o1.value_nz(p_map,key,cols,coeff);
-        o2.value_nz(p_map,key,cols,coeff);
+        o1.template value_nz<Sys_eqs>(p_map,key,cols,coeff,comp);
+        o2.template value_nz<Sys_eqs>(p_map,key,cols,coeff,comp);
     }
 
     /*! \brief Return the vector on which is acting
@@ -472,11 +472,11 @@ public:
 		return o1.value(key) / o2.value(key);
 	}
 
-    template<typename pmap_type, typename unordered_map_type, typename coeff_type>
-    inline void value_nz(pmap_type & p_map, const vect_dist_key_dx & key, unordered_map_type & cols, coeff_type & coeff) const
+    template<typename Sys_eqs,typename pmap_type, typename unordered_map_type, typename coeff_type>
+    inline void value_nz(pmap_type & p_map, const vect_dist_key_dx & key, unordered_map_type & cols, coeff_type & coeff, unsigned int comp) const
     {
-        o1.value_nz(p_map,key,cols,coeff);
-        o2.value_nz(p_map,key,cols,coeff);
+        o1.template value_nz<Sys_eqs>(p_map,key,cols,coeff,comp);
+        o2.template value_nz<Sys_eqs>(p_map,key,cols,coeff,comp);
     }
 
     /*! \brief Return the vector on which is acting
@@ -628,11 +628,11 @@ public:
 		return -(o1.value(key));
 	}
 
-    template<typename pmap_type, typename unordered_map_type, typename coeff_type>
-    inline void value_nz(pmap_type & p_map, const vect_dist_key_dx & key, unordered_map_type & cols, coeff_type & coeff) const
+    template<typename Sys_eqs, typename pmap_type, typename unordered_map_type, typename coeff_type>
+    inline void value_nz(pmap_type & p_map, const vect_dist_key_dx & key, unordered_map_type & cols, coeff_type & coeff, unsigned int comp) const
     {
 	    coeff_type coeff_tmp = -coeff;
-        o1.value_nz(p_map,key,cols,coeff_tmp);
+        o1.template value_nz<Sys_eqs>(p_map,key,cols,coeff_tmp, comp);
     }
 };
 
@@ -780,10 +780,10 @@ public:
 		return v;
 	}
 
-    template<typename pmap_type, typename unordered_map_type, typename coeff_type>
-    inline void value_nz(pmap_type & p_map, const vect_dist_key_dx & key, unordered_map_type & cols, coeff_type & coeff) const
+    template<typename Sys_eqs, typename pmap_type, typename unordered_map_type, typename coeff_type>
+    inline void value_nz(pmap_type & p_map, const vect_dist_key_dx & key, unordered_map_type & cols, coeff_type & coeff, unsigned int comp) const
     {
-	    cols[p_map. template getProp<0>(key)] += coeff;
+	    cols[p_map. template getProp<0>(key)*Sys_eqs::nvar + comp] += coeff;
     }
 
     inline vector_dist_expression_op<vector_dist_expression<prp,vector>,boost::mpl::int_<1>,VECT_COMP> operator[](int comp)
@@ -977,10 +977,34 @@ public:
 		return get_vector_dist_expression_op<n,n == rank_gen<property_act>::type::value>::get(o1,key,comp);
 	}
 
-    template<typename pmap_type, typename unordered_map_type, typename coeff_type>
-    inline void value_nz(pmap_type & p_map, const vect_dist_key_dx & key, unordered_map_type & cols, coeff_type & coeff) const
+	/*! \brief Return the result of the expression
+	 *
+	 * \note this function must be deactivated on transitional objects. Suppose we are slicing a tensor of rank 2
+	 *            an object of rank 1 is implicitly created for such object we have to deactivate this function
+	 *            because ill-formed
+	 *
+	 * \param key point where to evaluate
+	 *
+	 *
+	 */
+	inline auto get(const vect_dist_key_dx & key) const -> decltype(value(key))
+	{
+		return this->value(key);
+	}
+
+    template<typename Sys_eqs, typename pmap_type, typename unordered_map_type, typename coeff_type>
+    inline void value_nz(pmap_type & p_map, const vect_dist_key_dx & key, unordered_map_type & cols, coeff_type & coeff, unsigned int comp_) const
     {
-        o1.value_nz(p_map,key,cols,coeff,comp);
+#ifdef SE_CLASS1
+
+    	if (n != 1)
+    	{
+    		std::cout << __FILE__ << ":" << __LINE__ << " Error it only work for tensore of rank 1 ... like vectors " << std::endl;
+    	}
+
+#endif
+
+        o1.template value_nz<Sys_eqs>(p_map,key,cols,coeff,comp_ + comp[0]);
     }
 
     inline vector_dist_expression_op<exp1,boost::mpl::int_<2>,VECT_COMP> operator[](int comp_)
@@ -1130,10 +1154,12 @@ public:
 	{
 		return d;
 	}
-    template<typename pmap_type, typename unordered_map_type, typename coeff_type>
-    inline void value_nz(pmap_type & p_map, const vect_dist_key_dx & key, unordered_map_type & cols, coeff_type & coeff) const
+
+
+    template<typename Sys_eqs, typename pmap_type, typename unordered_map_type, typename coeff_type>
+    inline void value_nz(pmap_type & p_map, const vect_dist_key_dx & key, unordered_map_type & cols, coeff_type & coeff, unsigned int comp) const
     {
-        cols[p_map. template getProp<0>(key)] += coeff;
+        cols[p_map. template getProp<0>(key)*Sys_eqs::nvar + comp] += coeff;
     }
 };
 
-- 
GitLab