diff --git a/mpi_include b/mpi_include
index 16981421116cb2ce7f9dcf3dba32df467827856c..4439b1ca3ff6f7e0872a08610a3d8f99ee1db3a9 100644
--- a/mpi_include
+++ b/mpi_include
@@ -1 +1 @@
--I
\ No newline at end of file
+-I/home/i-bird/MPI/include
\ No newline at end of file
diff --git a/mpi_libs b/mpi_libs
index 0519ecba6ea913e21689ec692e81e9e4973fbf73..3a37a5cd2100d56757bcde4ef6ae508b3f7c0c68 100644
--- a/mpi_libs
+++ b/mpi_libs
@@ -1 +1 @@
-
\ No newline at end of file
+-Wl,-rpath -Wl,/home/i-bird/MPI/lib -Wl,--enable-new-dtags -pthread /home/i-bird/MPI/lib/libmpi.so
\ No newline at end of file
diff --git a/src/DCPSE_op/DCPSE_copy/Dcpse.hpp b/src/DCPSE_op/DCPSE_copy/Dcpse.hpp
index 875dcc4a8c69a893cfff505d8316d814470d65ab..d24f97a21dcee46b9916b82963d0d1bfd430846a 100644
--- a/src/DCPSE_op/DCPSE_copy/Dcpse.hpp
+++ b/src/DCPSE_op/DCPSE_copy/Dcpse.hpp
@@ -15,6 +15,26 @@
#include "DcpseDiagonalScalingMatrix.hpp"
#include "DcpseRhs.hpp"
+template<bool cond>
+struct is_scalar {
+ template<typename op_type>
+ static auto
+ analyze(const vect_dist_key_dx &key, op_type &o1) -> typename std::remove_reference<decltype(o1.value(
+ key))>::type {
+ return o1.value(key);
+ };
+};
+
+template<>
+struct is_scalar<false> {
+ template<typename op_type>
+ static auto
+ analyze(const vect_dist_key_dx &key, op_type &o1) -> typename std::remove_reference<decltype(o1.value(
+ key))>::type {
+ return o1.value(key);
+ };
+};
+
template<unsigned int dim, typename vector_type>
class Dcpse {
public:
@@ -219,26 +239,6 @@ public:
}
- template<bool cond>
- struct is_scalar {
- template<typename op_type>
- static auto
- analyze(const vect_dist_key_dx &key, op_type &o1) -> typename std::remove_reference<decltype(o1.value(
- key))>::type {
- return o1.value(key);
- };
- };
-
- template<>
- struct is_scalar<false> {
- template<typename op_type>
- static auto
- analyze(const vect_dist_key_dx &key, op_type &o1) -> typename std::remove_reference<decltype(o1.value(
- key))>::type {
- return o1.value(key);
- };
- };
-
/*! \brief Get the number of neighbours
*
* \return the number of neighbours
@@ -300,16 +300,13 @@ public:
}
/**
- * Computes the value of the differential operator on all the particles,
- * using the f values stored at the fValuePos position in the aggregate
- * and storing the resulting Df values at the DfValuePos position in the aggregate.
- * @tparam fValuePos Position in the aggregate of the f values to use.
- * @tparam DfValuePos Position in the aggregate of the Df values to store.
- * @param particles The set of particles to iterate over.
+ * Computes the value of the differential operator for one particle for o1 representing a scalar
+ *
+ * \param key particle
+ * \param o1 source property
+ * \return the selected derivative
+ *
*/
-
-
-
template<typename op_type>
auto computeDifferentialOperator(const vect_dist_key_dx &key,
op_type &o1) -> decltype(is_scalar<std::is_fundamental<decltype(o1.value(
@@ -317,8 +314,6 @@ public:
typedef decltype(is_scalar<std::is_fundamental<decltype(o1.value(key))>::value>::analyze(key, o1)) expr_type;
- //typedef typename decltype(o1.value(key))::blabla blabla;
-
T sign = 1.0;
if (differentialOrder % 2 == 0) {
sign = -1;
@@ -347,6 +342,15 @@ public:
return Dfxp;
}
+ /**
+ * Computes the value of the differential operator for one particle for o1 representing a vector
+ *
+ * \param key particle
+ * \param o1 source property
+ * \param i component
+ * \return the selected derivative
+ *
+ */
template<typename op_type>
auto computeDifferentialOperator(const vect_dist_key_dx &key,
op_type &o1,
diff --git a/src/DCPSE_op/DCPSE_copy/DcpseRhs.hpp b/src/DCPSE_op/DCPSE_copy/DcpseRhs.hpp
index c11e0dcd27ddbe15d4fe8d2b1e55e6aee03c30ed..b3ff298778bdb6194f589dfb900f29309b99e285 100644
--- a/src/DCPSE_op/DCPSE_copy/DcpseRhs.hpp
+++ b/src/DCPSE_op/DCPSE_copy/DcpseRhs.hpp
@@ -53,7 +53,7 @@ MatrixType &DcpseRhs<dim>::getVector(MatrixType &b)
if (b(0,0) == 0.0 && sign == 1)
{
- b(0,0) = 0.0;
+ b(0,0) = 10.0;
}
return b;
diff --git a/src/DCPSE_op/DCPSE_op_test.cpp b/src/DCPSE_op/DCPSE_op_test.cpp
index 0d9d5252af37f80633634cb843050dd3ace55ca1..c84d6eaf29b4a925c9d69e84e00363540bbf8f65 100644
--- a/src/DCPSE_op/DCPSE_op_test.cpp
+++ b/src/DCPSE_op/DCPSE_op_test.cpp
@@ -605,7 +605,7 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests)
if (up.isInside(xp) == true) {
up_p.add();
up_p.last().get<0>() = p.getKey();
- Particles.getProp<1>(p) = -2/sz[0];
+ Particles.getProp<1>(p) = -2.0/sz[0];
}
else if (down.isInside(xp) == true) {
dw_p.add();
@@ -625,11 +625,11 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests)
}
Particles.write_frame("Re1000-1e-3-Lid_sf",0);
- Derivative_x Dx(Particles, 2, rCut,3);
- Derivative_y Dy(Particles, 2, rCut,3);
- Gradient Grad(Particles, 2, rCut,3);
- Laplacian Lap(Particles, 2, rCut, 3);
- Curl2D Curl(Particles, 2, rCut, 3);
+ Derivative_x Dx(Particles, 1, rCut,3);
+ Derivative_y Dy(Particles, 1, rCut,3);
+ Gradient Grad(Particles, 1, rCut,3);
+ Laplacian Lap(Particles, 1, rCut, 3);
+ Curl2D Curl(Particles, 1, rCut, 3);
DCPSE_scheme<equations,decltype(Particles)> Solver(2*rCut, Particles);
auto Sf_Poisson = Lap(Sf);
Solver.impose(Sf_Poisson, bulk, prop_id<1>());
@@ -639,9 +639,13 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests)
Solver.impose(Sf, r_p, 0);
Solver.solve(Sf);
V=Curl(Sf);
+
+ Particles.write("Initialization");
+
+
double dt=0.003;
double nu=0.01;
- int n=2;
+ int n=150;
std::cout<<"Init Done"<<std::endl;
for(int i =1; i<=n ;i++)
{ dW=nu*Lap(W)+Dx(Sf)*Dy(W)-Dy(Sf)*Dx(W);
@@ -1421,6 +1425,71 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests)
}
+
+ BOOST_AUTO_TEST_CASE(dcpse_slice) {
+ const size_t sz[2] = {31,31};
+ Box<2, double> box({0, 0}, {1,1});
+ size_t bc[2] = {NON_PERIODIC, NON_PERIODIC};
+ double spacing = box.getHigh(0) / (sz[0] - 1);
+ Ghost<2, double> ghost(spacing * 3);
+ double rCut = 2.0 * spacing;
+
+ vector_dist<2, double, aggregate<double,VectorS<2, double>,double>> Particles(0, box, bc, ghost);
+
+ auto it = Particles.getGridIterator(sz);
+ while (it.isNext()) {
+ Particles.add();
+ auto key = it.get();
+ double x = key.get(0) * it.getSpacing(0);
+ Particles.getLastPos()[0] = x;
+ double y = key.get(1) * it.getSpacing(1);
+ Particles.getLastPos()[1] = y;
+
+ Particles.getLastProp<1>()[0] = sin(x+y);
+ Particles.getLastProp<1>()[1] = cos(x+y);
+
+ ++it;
+ }
+
+ Particles.map();
+ Particles.ghost_get<0>();
+
+
+ auto P = getV<0>(Particles);
+ auto V = getV<1>(Particles);
+ auto S = getV<2>(Particles);
+
+ Derivative_x Dx(Particles, 2, rCut,2);
+
+ P = Dx(V[0]);
+ S = V[0]*V[0] + V[1]*V[1];
+
+ auto it2 = Particles.getDomainIterator();
+
+ double err = 0.0;
+
+ while (it2.isNext())
+ {
+ auto p = it2.get();
+
+ if (fabs(Particles.getProp<0>(p) - Particles.getProp<1>(p)[1]) >= err )
+ {
+ err = fabs(Particles.getProp<0>(p) - Particles.getProp<1>(p)[1]);
+ }
+
+ if (fabs(Particles.getProp<2>(p) - 1.0) >= err )
+ {
+ err = fabs(Particles.getProp<2>(p) - 1.0);
+ }
+
+ ++it2;
+ }
+
+ std::cout << "ERR " << err << std::endl;
+
+ Particles.write("test_out");
+ }
+
BOOST_AUTO_TEST_SUITE_END()
diff --git a/src/Operators/Vector/vector_dist_operators.hpp b/src/Operators/Vector/vector_dist_operators.hpp
index df7001d6acbb44e7fb3016dcaf93c54b7e1397ca..62ae5170fc2deb6950fb105d520d4081521da764 100644
--- a/src/Operators/Vector/vector_dist_operators.hpp
+++ b/src/Operators/Vector/vector_dist_operators.hpp
@@ -62,16 +62,20 @@
#define VECT_ROUND 88
#define VECT_NEARBYINT 89
#define VECT_RINT 90
+#define VECT_PMUL 91
+#define VECT_SUB_UNI 92
+#define VECT_SUM_REDUCE 93
+#define VECT_COMP 94
+
+
#define VECT_DCPSE 100
#define VECT_DCPSE_V 101
#define VECT_DCPSE_V_SUM 102
#define VECT_DCPSE_V_DOT 103
#define VECT_DCPSE_V_DIV 104
#define VECT_DCPSE_V_CURL2D 105
-#define VECT_PMUL 91
-#define VECT_SUB_UNI 92
-#define VECT_SUM_REDUCE 93
+
template<bool cond, typename exp1, typename exp2>
struct first_or_second
@@ -587,6 +591,8 @@ class vector_dist_expression_op<exp1,void,VECT_SUB_UNI>
public:
+ typedef typename exp1::vtype vtype;
+
//! constructor from an expresssion
vector_dist_expression_op(const exp1 & o1)
:o1(o1)
@@ -598,6 +604,19 @@ public:
o1.init();
}
+ /*! \brief Return the vector on which is acting
+ *
+ * It return the vector used in getVExpr, to get this object
+ *
+ * \return the vector
+ *
+ */
+ const vtype & getVector() const
+ {
+ return o1.getVector();
+ }
+
+
//! return the result of the expression
template<typename r_type=typename std::remove_reference<decltype(-(o1.value(vect_dist_key_dx(0))))>::type > inline r_type value(const vect_dist_key_dx & key) const
{
@@ -612,6 +631,59 @@ public:
}
};
+/*! \brief it take an expression and create the negatove of this expression
+ *
+ *
+ */
+template <typename exp1>
+class vector_dist_expression_op<exp1,void,VECT_COMP>
+{
+ //! expression 1
+ const exp1 o1;
+
+ //! component
+ int comp;
+
+public:
+
+ typedef typename exp1::vtype vtype;
+
+ //! constructor from an expresssion
+ vector_dist_expression_op(const exp1 & o1, int comp)
+ :o1(o1),comp(comp)
+ {}
+
+ /*! \brief Return the vector on which is acting
+ *
+ * It return the vector used in getVExpr, to get this object
+ *
+ * \return the vector
+ *
+ */
+ const vtype & getVector() const
+ {
+ return o1.getVector();
+ }
+
+ //! initialize the expression tree
+ inline void init() const
+ {
+ o1.init();
+ }
+
+ //! return the result of the expression
+ template<typename r_type=typename std::remove_reference<decltype(o1.value(vect_dist_key_dx(0))[comp])>::type > inline r_type value(const vect_dist_key_dx & key) const
+ {
+ return o1.value(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
+ {
+ o1.value_nz(p_map,key,cols,coeff,comp);
+ }
+};
+
/*! \brief Main class that encapsulate a vector properties operand to be used for expressions construction
*
* \tparam prp property involved
@@ -759,6 +831,13 @@ public:
{
cols[p_map. template getProp<0>(key)] += coeff;
}
+
+ inline vector_dist_expression_op<vector_dist_expression<prp,vector>,void,VECT_COMP> operator[](int comp)
+ {
+ vector_dist_expression_op<vector_dist_expression<prp,vector>,void,VECT_COMP> v_exp(*this,comp);
+
+ return v_exp;
+ }
};
/*! \Create an expression from a vector property