main.cpp

#include "Vector/vector_dist.hpp"
#include "Operators/Vector/vector_dist_operators.hpp"
 
 
constexpr int A = 0;
constexpr int B = 1;
constexpr int C = 2;
constexpr int D = 3;
 
 
 
 
// Kernel are structures
struct exp_kernel
{
    // variance of the exponential
    float sigma;
 
    // Constructor require to define the variance
    exp_kernel(float sigma)
    :sigma(sigma)
    {}
 
    // The kernel function itself K(x_p,x_q,A_p,A_q). The name MUST be value and require 4 arguments.
    // p,q,A_p,A_q Position of p, position of q, value of the property A on p, value of the property A on q
    // and it return the same type of the property A
    inline Point<3,double> value(const Point<3,float> & p, const Point<3,float> & q,Point<3,double> & Pp,Point<3,double> & Pq)
    {
        // calculate the distance between p and q
        float dist = norm(p-q);
 
        // Calculate the exponential
        return Pq * exp(- dist * dist / sigma);
    }
 
    // The kernel function itself K(p,q,A_p,A_q,particles). The name MUST be value and require 5 arguments.
    // p,q,A_p,A_q,particles size_t index of p, size_t index of q, value of the property A on p, value of the property A on q,
    // original vector of the particles p, and it return the same type of the property A
    inline Point<3,double> value(size_t p, size_t q, Point<3,double> & Pp, Point<3,double> & Pq, const vector_dist<3,double, aggregate<double,double,Point<3,double>,Point<3,double>> > & v)
    {
        // calculate the distance between p and q
        float dist = norm(p-q);
 
        // Calculate the exponential
        return Pq * exp(- dist * dist / sigma);
    }
};
 
 
int main(int argc, char* argv[])
{
 
    // initialize the library
    openfpm_init(&argc,&argv);
 
    // Here we define our domain a 2D box with intervals from 0 to 1.0
    Box<3,double> domain({0.0,0.0,0.0},{1.0,1.0,1.0});
 
    // Here we define the boundary conditions of our problem
    size_t bc[3]={PERIODIC,PERIODIC,PERIODIC};
 
    // extended boundary around the domain, and the processor domain
    Ghost<3,double> g(0.01);
 
    // Delta t
    double dt = 0.01;
 
 
 
    // distributed vectors
    vector_dist<3,double, aggregate<double,double,Point<3,double>,Point<3,double>> > vd(4096,domain,bc,g);
    vector_dist<3,double, aggregate<double> > vd2(4096,domain,bc,g);
 
 
 
    // Assign random position to the vector dist
    auto it = vd.getDomainIterator();
 
    while (it.isNext())
    {
        auto p = it.get();
 
        vd.getPos(p)[0]  = (double)rand() / RAND_MAX;
        vd.getPos(p)[1]  = (double)rand() / RAND_MAX;
        vd.getPos(p)[2]  = (double)rand() / RAND_MAX;
 
        vd2.getPos(p)[0] = (double)rand() / RAND_MAX;
        vd2.getPos(p)[1]  = (double)rand() / RAND_MAX;
        vd2.getPos(p)[2]  = (double)rand() / RAND_MAX;
 
        ++it;
    }
 
 
 
    vd.map();
    vd.map();
 
 
 
    // vA is an alias for the property A of the vector vd
    // vB is an alias for the property V of the vector vd
    // vC is ...
    auto vA = getV<A>(vd);
    auto vB = getV<B>(vd);
    auto vC = getV<C>(vd);
    auto vD = getV<D>(vd);
 
    // This indicate the particle position
    auto vPOS = getV<PROP_POS>(vd);
 
    // same concept for the vector v2
    auto v2A = getV<0>(vd);
    auto v2POS = getV<PROP_POS>(vd);
 
 
 
    // Assign 1 to the property A
    vA = 1;
 
    // Equal the property A to B
    vB = vA;
 
    // Assign to the property C and for each component 4.0
    vC = 4;
 
    // Assign the position of the particle to the property D
    vD = vPOS;
 
 
 
    // In vA we set the distance from the origin for each particle
    vA = norm(vPOS);
 
    //
    // For each particle p calculate the expression under
    //
    // NOTE sin(2.0 * vD)  and exp(5.0 * vD) are both vector
    //
    //              so this is a scalar product
    //                      |
    //                      V
    vA = vA + sin(2.0 * vD) * exp(1.2 * vD);
    //       |_____________________________|
    //                      |
    //                      V
    //              and this is a number
    //
    //
 
 
    // pmul indicate component-wise multiplication the same as .* in Matlab
    vC = pmul(vD,vD) * dt;
 
    // Normalization of the vector vD
    vD = vD / sqrt( vD * vD );
 
 
 
    // we write the result on file
    vd.write("output");
 
 
 
    // Constant point
    Point<3,double> p0({0.5,0.5,0.5});
 
    // we cannot use p0 directly we have to create an expression from it
    auto p0_e = getVExpr(p0);
 
    // here we are creating an expression, expr1 collect the full expression but does not produce any
    // code execution
    auto expr1 = 1.0/2.0/sqrt(M_PI)*exp(-(v2POS-p0_e)*(v2POS-p0_e)/2.0);
 
    // here the expression is executed on all particles of v2 and the result assigned to the property A of v2
    v2A = expr1;
 
 
 
    assign(vA,1.0,               // vA = 1
           vB,expr1,          // vB = 1.0/2.0/sqrt(M_PI)*exp(-(v2POS-p0_e)*(v2POS-p0_e)/2.0)
           vC,vD/norm(vD),    // vC = vD/norm(vD)
           vD,2.0*vC);        // vD = 2.0*vC // here vC = vD/norm(vD)
 
 
 
    // Create a kernel with sigma 0.5
    exp_kernel ker(0.5);
 
    // Get a Cell list with 0.1 of cutoff-radius
    auto cl = vd.getCellList(0.1);
 
    // We are going to do some calculation that require ghost, so sync it 3 == vD
    vd.ghost_get<3>();
 
 
 
    vC = applyKernel_in( vD / norm(vD) ,vd,cl,ker) + vD;
 
 
 
    // Second version in this case the kernel is called with the particles id for p
    // and
    vC = applyKernel_in_gen( vD / norm(vD) ,vd,cl,ker) + vD;
 
 
    // ok, because does not use applyKernel
    vD = vD / norm(vD);
 
 
 
    vd.deleteGhost();
    vd.write("output2");
 
 
 
    openfpm_finalize();
 
}