diff --git a/src/DCPSE_op/DCPSE_Solver.hpp b/src/DCPSE_op/DCPSE_Solver.hpp
index 6aaa290943f891640ad193b35fedd6fa67a96a74..e8ee6419082c5b5f5d532a5e87108f61c43cda6f 100644
--- a/src/DCPSE_op/DCPSE_Solver.hpp
+++ b/src/DCPSE_op/DCPSE_Solver.hpp
@@ -278,47 +278,60 @@ class DCPSE_scheme: public MatrixAssembler
}
}
-public:
-
-/* void * dcpse_solver;
-
-public:
-
- template<typename particles_type>
- Solver(expr type,particles_type & parts, unsigned int ord ,typename particles_type::stype rCut,double scaling_factor=dcpse_scaling_factor)
+ /*! \brief Solve an equation
+ *
+ * \warning exp must be a scalar type
+ *
+ * \param exp where to store the result
+ *
+ */
+ template<typename solType, typename expr_type>
+ void copy_impl(solType & x, expr_type exp, unsigned int comp)
{
- typedef Dcpse<particles_type::dims,particles_type> DCPSE_type;
-
- dcpse_solver = new unsigned char [particles_type::dims*sizeof(DCPSE_type)];
-
- SparseMatrix<double, int, PETSC_BASE> MAT(sz[0]*sz[1], sz[0]*sz[1], sz[0]*sz[1]);
+ auto parts = exp.getVector();
- Dcpse<particles_type::dims,particles_type> * dcpse_ptr = (Dcpse<particles_type::dims,particles_type> *)dcpse;
+ auto it = parts.getDomainIterator();
- for (int i = 0 ; i < particles_type::dims ; i++)
+ while (it.isNext())
{
-
+ auto p = it.get();
+ exp.value(p) = x(p.getKey()*Sys_eqs::nvar + comp);
+ ++it;
}
}
- template<typename operand_type>
- vector_dist_expression_op<operand_type,Dcpse<operand_type::vtype::dims,typename operand_type::vtype>,VECT_DCPSE_V> operator()(operand_type arg)
+ template<typename solType, typename exp1, typename ... othersExp>
+ void copy_nested(solType & x, unsigned int & comp, exp1 exp, othersExp ... exps)
{
- typedef Dcpse<operand_type::vtype::dims,typename operand_type::vtype> dcpse_type;
+ copy_impl(x,exp,comp);
+ comp++;
+
+ copy_nested(x,comp,exps ...);
+ }
- return vector_dist_expression_op<operand_type,dcpse_type,VECT_DCPSE_V>(arg,*(dcpse_type(*)[operand_type::vtype::dims])dcpse);
- }*/
+ template<typename solType, typename exp1>
+ void copy_nested(solType & x, unsigned int & comp, exp1 exp)
+ {
+ copy_impl(x,exp,comp);
+ comp++;
+ }
+
+public:
+
+ /*! \brief Solve an equation
+ *
+ * \warning exp must be a scalar type
+ *
+ * \param exp where to store the result
+ *
+ */
template<typename expr_type>
void solve(expr_type exp)
{
umfpack_solver<double> solver;
auto x = solver.solve(getA(opt),getB(opt));
- //petsc_solver<double> solver;
- //solver.setSolver(KSPBCGS);
- //solver.setMaxIter(1000);
- //solver.setRestart(200);
- //auto x = solver.try_solve(getA(),getB());
+
auto parts = exp.getVector();
auto it = parts.getDomainIterator();
@@ -331,6 +344,30 @@ public:
}
}
+ /*! \brief Solve an equation
+ *
+ * \warning exp must be a scalar type
+ *
+ * \param exp where to store the result
+ *
+ */
+ template<typename ... expr_type>
+ void solve2(expr_type ... exps)
+ {
+ std::cout << sizeof...(exps) << " " << Sys_eqs::nvar << std::endl;
+
+ if (sizeof...(exps) != Sys_eqs::nvar)
+ {std::cerr << __FILE__ << ":" << __LINE__ << " Error the number of properties you gave does not match the solution in\
+ dimensionality, I am expecting " << Sys_eqs::nvar <<
+ " properties " << std::endl;};
+
+ umfpack_solver<double> solver;
+ auto x = solver.solve(getA(opt),getB(opt));
+
+ unsigned int comp = 0;
+ copy_nested(x,comp,exps ...);
+ }
+
/*! \brief Constructor
*
* The periodicity is given by the grid b_g
diff --git a/src/DCPSE_op/DCPSE_op_test.cpp b/src/DCPSE_op/DCPSE_op_test.cpp
index 66ea2e27e870f0000a20820334a0c5fe380487af..c4087ac68576e967de3e54a1eed4d22181a2c3dd 100644
--- a/src/DCPSE_op/DCPSE_op_test.cpp
+++ b/src/DCPSE_op/DCPSE_op_test.cpp
@@ -2445,7 +2445,7 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests)
Solver.impose(v[1], dw_p, RHS[1],vy);
Solver.impose(v[0], l_p, RHS[0],vx);
Solver.impose(v[1], l_p, RHS[1],vy);
- Solver.solve(sol);
+ Solver.solve2(sol[0],sol[1]);
DCPSE_sol=Lap(sol);
double worst1 = 0.0;
double worst2 = 0.0;
diff --git a/src/Operators/Vector/vector_dist_operators.hpp b/src/Operators/Vector/vector_dist_operators.hpp
index 64ee1efc9eaf3101cd09bb953fd11a76e50a20eb..b981d6fe65945874504153feb43002b61ad9cb21 100644
--- a/src/Operators/Vector/vector_dist_operators.hpp
+++ b/src/Operators/Vector/vector_dist_operators.hpp
@@ -802,7 +802,7 @@ template<unsigned int, bool is_valid>
struct get_vector_dist_expression_op
{
template<typename exp_type>
- inline static auto get(exp_type & o1, const vect_dist_key_dx & key) -> typename std::remove_reference<decltype(o1.value(vect_dist_key_dx(0)))>::type
+ inline static auto get(exp_type & o1, const vect_dist_key_dx & key) -> decltype(o1.value(vect_dist_key_dx(0)))
{
return o1.value(key);
}
@@ -844,7 +844,7 @@ template<>
struct get_vector_dist_expression_op<1,true>
{
template<typename exp_type>
- static auto get(exp_type & o1, const vect_dist_key_dx & key, const int (& comp)[1]) -> typename std::remove_reference<decltype(o1.value(vect_dist_key_dx(0))[0])>::type
+ static auto get(exp_type & o1, const vect_dist_key_dx & key, const int (& comp)[1]) -> decltype(o1.value(vect_dist_key_dx(0))[0])
{
return o1.value(key)[comp[0]];
}
@@ -866,7 +866,7 @@ template<>
struct get_vector_dist_expression_op<2,true>
{
template<typename exp_type>
- static auto get(exp_type & o1, const vect_dist_key_dx & key, const int (& comp)[2]) -> typename std::remove_reference<decltype(o1.value(vect_dist_key_dx(0))[0][0])>::type
+ static auto get(exp_type & o1, const vect_dist_key_dx & key, const int (& comp)[2]) -> decltype(o1.value(vect_dist_key_dx(0))[0][0])
{
return o1.value(key)[comp[0]][comp[1]];
}