src/DCPSE/DcpseDiagonalScalingMatrix.hpp

+//
+// Created by tommaso on 29/03/19.
+// Modified by Serhii
+//
+
+#ifndef OPENFPM_PDATA_DCPSEDIAGONALSCALINGMATRIX_HPP
+#define OPENFPM_PDATA_DCPSEDIAGONALSCALINGMATRIX_HPP
+
+#include "MonomialBasis.hpp"
+#include "Support.hpp"
+
+
+template <unsigned int dim, typename monomialBasis_type = MonomialBasis<dim>>
+class DcpseDiagonalScalingMatrix
+{
+private:
+    const monomialBasis_type& monomialBasis;
+
+public:
+    DcpseDiagonalScalingMatrix(const monomialBasis_type &monomialBasis) : monomialBasis(monomialBasis) {}
+
+    template <typename T, typename MatrixType, typename vector_type, typename vector_type2>
+    void buildMatrix(MatrixType &M, Support support, T eps, vector_type & particlesFrom , vector_type2 & particlesTo)
+    {
+        // Check that all the dimension constraints are met
+        assert(support.size() >= monomialBasis.size());
+        assert(M.rows() == support.size());
+        assert(M.cols() == support.size());
+
+        Point<dim,typename vector_type::stype> ref_p = particlesTo.getPosOrig(support.getReferencePointKey());
+
+        // Fill the diagonal matrix
+        M.setZero(); // Make sure the rest of the matrix is zero!
+        const auto& support_keys = support.getKeys();
+        size_t N = support_keys.size();
+        for (size_t i = 0; i < N; ++i)
+        {
+            const auto& pt = support_keys.get(i);
+        	Point<dim,typename vector_type::stype> p = ref_p;
+        	p -= particlesFrom.getPosOrig(pt);
+
+            M(i,i) = exp(- norm2(p) / (2.0 * eps * eps));
+        }
+    }
+
+    template <typename T, typename vector_type, typename vector_type2>
+    __host__ __device__ void buildMatrix(T* M, size_t supportRefKey, size_t supportKeysSize, const size_t* supportKeys, T eps, vector_type & particlesFrom, vector_type2 & particlesTo)
+    {
+        // Check that all the dimension constraints are met
+        assert(supportKeysSize >= monomialBasis.size());
+
+        Point<dim,typename vector_type::stype> ref_p = particlesTo.getPos(supportRefKey);
+
+        for (size_t i = 0; i < supportKeysSize; ++i)
+        {
+            size_t pt = supportKeys[i];
+            Point<dim,typename vector_type::stype> p = ref_p;
+            p -= particlesFrom.getPos(pt);
+
+            M[i] = exp(- norm2(p) / (2.0 * eps * eps));
+        }
+    }
+};
+
+#endif //OPENFPM_PDATA_DCPSEDIAGONALSCALINGMATRIX_HPP

src/DCPSE/DcpseInterpolation.hpp

+//
+// Created by Abhinav Singh on 03.11.21.
+//
+
+#ifndef OPENFPM_PDATA_DCPSEINTERPOLATION_HPP
+#define OPENFPM_PDATA_DCPSEINTERPOLATION_HPP
+#include "DCPSE/Dcpse.hpp"
+
+/*! \brief Class for Creating the DCPSE Operator For the function approximation objects and computes DCPSE Kernels.
+ *
+ *
+ * \param parts particle set
+ * \param ord order of convergence of the operator
+ * \param rCut Argument for cell list construction
+ * \param oversampling_factor multiplier to the minimum no. of particles required by the operator in support
+ * \param support_options default:N_particles, Radius can be used to select all particles inside rCut. Overrides oversampling.
+ *
+ * \return Operator Dx which is a function on Vector_dist_Expressions
+ *
+ */
+template<typename particlesFrom_type, typename particlesTo_type>
+class PPInterpolation 
+{
+
+    void *dcpse;
+
+    particlesFrom_type & particlesFrom;
+    particlesTo_type & particlesTo;
+
+public:
+    /*! \brief Constructor for Creating the DCPSE Operator Dx and objects and computes DCPSE Kernels.
+     *
+     *
+     * \param parts particle set
+     * \param ord order of convergence of the operator
+     * \param rCut Argument for cell list construction
+     * \param oversampling_factor multiplier to the minimum no. of particles required by the operator in support
+     * \param support_options default:N_particles, Radius can be used to select all particles inside rCut. Overrides oversampling.
+     *
+     * \return Operator F which is a function on Vector_dist_Expressions
+     *
+     */
+    PPInterpolation(particlesFrom_type &particlesFrom,particlesTo_type &particlesTo, unsigned int ord, typename particlesFrom_type::stype rCut,
+                      double oversampling_factor = dcpse_oversampling_factor,
+                      support_options opt = support_options::RADIUS)
+    :particlesFrom(particlesFrom),particlesTo(particlesTo)
+    {
+        Point<particlesFrom_type::dims, unsigned int> p;
+        p.zero();
+        dcpse = new Dcpse<particlesFrom_type::dims, particlesFrom_type,particlesTo_type>(particlesFrom,particlesTo, p, ord, rCut, oversampling_factor, opt);
+    }
+
+    void deallocate() {
+        delete (Dcpse<particlesFrom_type::dims, particlesFrom_type, particlesTo_type> *) dcpse;
+    }
+
+   /* template<typename operand_type>
+    vector_dist_expression_op<operand_type, Dcpse_type<operand_type::vtype::dims, typename operand_type::vtype>, VECT_DCPSE>
+    operator()(operand_type arg) {
+        typedef Dcpse_type<operand_type::vtype::dims, typename operand_type::vtype> dcpse_type;
+        return vector_dist_expression_op<operand_type, dcpse_type, VECT_DCPSE>(arg, *(dcpse_type *) dcpse);
+    }*/
+
+   template<unsigned int prp1,unsigned int prp2>
+   void p2p() {
+       auto dcpse_temp = (Dcpse<particlesFrom_type::dims, particlesFrom_type, particlesTo_type>*) dcpse;
+       dcpse_temp->template p2p<prp1,prp2>();
+
+   }
+
+    // template<unsigned int prp, typename particles_type>
+    // void DrawKernel(particles_type &particles, int k) {
+    //     auto dcpse_temp = (Dcpse_type<particlesFrom_type::dims, particlesFrom_type, particlesTo_type> *) dcpse;
+    //     dcpse_temp->template DrawKernel<prp>(particles, k);
+
+    // }
+
+    // template<unsigned int prp, typename particles_type>
+    // void DrawKernelNN(particles_type &particles, int k) {
+    //     auto dcpse_temp = (Dcpse_type<particlesFrom_type::dims, particlesFrom_type,particlesTo_type> *) dcpse;
+    //     dcpse_temp->template DrawKernelNN<prp>(particles, k);
+
+    // }
+
+    // template<typename particles_type>
+    // void checkMomenta(particles_type &particles) {
+    //     auto dcpse_temp = (Dcpse_type<particles_type::dims, particlesFrom_type, particlesTo_type> *) dcpse;
+    //     dcpse_temp->checkMomenta(particles);
+
+    // }
+
+    /*! \brief Method for Updating the DCPSE Operator by recomputing DCPSE Kernels.
+     *
+     *
+     * \param parts particle set
+     */
+    void update() {
+        auto dcpse_temp = (Dcpse<particlesFrom_type::dims, particlesFrom_type, particlesTo_type> *) dcpse;
+        dcpse_temp->initializeUpdate(particlesFrom,particlesTo);
+
+    }
+
+};
+
+
+
+#endif //OPENFPM_PDATA_DCPSEINTERPOLATION_HPP

src/DCPSE/DcpseRhs.hpp

+//
+// Created by tommaso on 28/03/19.
+//
+
+#ifndef OPENFPM_PDATA_DCPSERHS_HPP
+#define OPENFPM_PDATA_DCPSERHS_HPP
+
+#include "MonomialBasis.hpp"
+
+template<unsigned int dim>
+class DcpseRhs
+{
+private:
+    const Point<dim, unsigned int> differentialSignature;
+    const MonomialBasis<dim> derivatives;
+    int sign;
+public:
+    DcpseRhs(const MonomialBasis<dim> &monomialBasis, const Point<dim, unsigned int> &differentialSignature);
+
+    template<typename T, typename MatrixType>
+    MatrixType &getVector(MatrixType &b);
+};
+
+// Definitions below
+
+template<unsigned int dim>
+DcpseRhs<dim>::DcpseRhs(const MonomialBasis<dim> &monomialBasis,
+                        const Point<dim, unsigned int> &differentialSignature)
+        : differentialSignature(differentialSignature), derivatives(monomialBasis.getDerivative(differentialSignature))
+{
+    unsigned int order = (Monomial<dim>(differentialSignature)).order();
+    if (order % 2 == 0)
+    {
+        sign = 1;
+    } else
+    {
+        sign = -1;
+    }
+}
+
+template<unsigned int dim>
+template<typename T, typename MatrixType>
+MatrixType &DcpseRhs<dim>::getVector(MatrixType &b)
+{
+    // The given vector V should have the right dimensions
+    assert(b.cols() == 1);
+    assert(b.rows() == derivatives.size());
+    for (unsigned int i = 0; i < derivatives.size(); ++i)
+    {
+        const Monomial<dim> dm = derivatives.getElement(i);
+        b(i, 0) = sign * dm.evaluate(Point<dim, T>(0));
+    }
+    //Choosing a(0,0) for even order 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) = 10.0;
+    }*/
+
+    return b;
+}
+
+#endif //OPENFPM_PDATA_DCPSERHS_HPP

src/Draw/DrawParticles_unit_tests.hpp

@@ -30,7 +30,7 @@ BOOST_AUTO_TEST_CASE(point_iterator)
 	size_t bc[3]={NON_PERIODIC,NON_PERIODIC,NON_PERIODIC};
 
 	// ghost, big enough to contain the interaction radius
-	Ghost<3,float> ghost(0.01);
+	Ghost<3,double> ghost(0.01);
 
 	vector_dist<3,double,aggregate<double>> vd(0,domain,bc,ghost);
 
@@ -87,7 +87,7 @@ BOOST_AUTO_TEST_CASE(point_iterator_skin)
 	size_t bc[3]={NON_PERIODIC,NON_PERIODIC,NON_PERIODIC};
 
 	// ghost, big enough to contain the interaction radius
-	Ghost<3,float> ghost(0.01);
+	Ghost<3,double> ghost(0.01);
 
 	vector_dist<3,double,aggregate<double>> vd(0,domain,bc,ghost);