1
votes

I have a particle system and I am using boost::geometry to approximate my elliptical particles as polygons and then use the intersection function of the library to find the overlap area. I am calculating an "inner" and "outer" ellipse(polygon) area to assign a "potential" for each particle-particle interaction.

My potential function is this:

    double Potential(Cell* current, Cell* next)
{
    double areaRep, areaAtt;
    double distance = Distance(current,next);
    double A1 = current->getLength();
    double B1 = A1/2.0;
    double theta1 = current->getTheta(); //*180.0/M_PI
    double x1 = current->getCurrX();
    double y1 = current->getCurrY();
    double A2 = next->getLength();
    double B2 = A2/2.0;
    double theta2 = next->getTheta();
    double x2 = next->getCurrX();
    double y2 = next->getCurrY();
    polygon_2d poly1, poly2, poly3, poly4;
    double lamda1, lamda2;
    lamda1 = 0.0005; lamda2 = 0.00001;
    if(distance < 2.0*1.5*A1) {
        ellipse2poly(theta1, A1, B1, x1, y1, &poly1);
        ellipse2poly(theta2, A2, B2, x2, y2, &poly2);
        areaRep = getOverlapingAreaPoly(poly1,poly2);
        ellipse2poly(theta1, 1.5*A1, 1.5*B1, x1, y1, &poly3);
        ellipse2poly(theta2, 1.5*A2, 1.5*B2, x2, y2, &poly4);
        areaAtt = getOverlapingAreaPoly(poly3, poly4);
        return (lamda1*areaRep - lamda2*areaAtt); 
    }
    else
        return 0.0;
}

The "polygonizing" function is:

int ellipse2poly(double theta, double A1, double B1, double H1, double K1, polygon_2d *po)
{
    using namespace boost::geometry;
    polygon_2d  poly;
    const int n = 20;
    double angle = theta; // cell orientation
    double a = A1; // Long semi-axis length
    double b = B1; // short semi-axis length
    double xc = H1; // current X position
    double yc = K1; // current Y position

    if(!n)
    {
        std::cout << "error ellipse(): n should be >0\n" <<std::endl;
        return 0;
    }
    double t = 0;
    int i = 0;
    double coor[2*n+1][2];

    double x, y;
    double step = M_PI/(double)n;
    double sinphi = sin(angle);
    double cosphi = cos(angle);
    for(i=0; i<2*n+1; i++)
    {   
        x = xc + a*cos(t)*cosphi - b*sin(t)*sinphi;
        y = yc + a*cos(t)*sinphi + b*sin(t)*cosphi;
        coor[i][0] = x;
        coor[i][1] = y;
        t += step;
    }
    assign_points(poly, coor);
    correct(poly);
    *po = poly;
    return 1;
}

And the returned area is:

double getOverlapingAreaPoly(polygon_2d poly, polygon_2d poly2)
{
    point_2d cent; //centre of overlaping area
    double overAreaPoly = 0.0;
    typedef std::vector<polygon_2d > polygon_list;
    polygon_list v;

    intersection(poly,poly2,v);
    for (polygon_list::const_iterator it = v.begin(); it != v.end(); ++it)
    {
        centroid(*it, cent);
        overAreaPoly = area(*it);
    }

    return overAreaPoly;
}

The function is called for every cell (particle) as long as it is not for the same one. Previously, using another method, one iteration of my algorithm would take approximately 43 ms for one iteration for 100 particles. Now it takes approximately 1 min(!!!), so I guess I have done something horribly wrong!

I have tested this only in MSVC2012 under win7 64bit. I will report back for Linux Mint with Qt 4.7.4.

EDIT: I have tested on Linux Mint with Qt 4.7.4 and it is running very reasonably; maybe 90-100 ms per iteration which is fine. I don't know what is wrong in win7...

1
Is your Potential() function fixed or are you willing to modify it? Can you provide an example for how you intend to use the potential?Matthias
Actually I fixed it. I just started a new project in Visual Studio, copied the source and header files, recompiled and it works very reasonably again! I don't know why that is actually... To answer your question, I am using a normalized vector of the distance between the particles to determine the direction of the interaction and then I multiply with the potential.Dima1982
Thank you for the clarification. Good to hear that it does work fine now :-)Matthias

1 Answers

1
votes

I have actually fixed it. I started a new project in Visual Studio and copied all source and header files, recompiled and everything runs smoothly now. I guess radically changing code and adding / subtracting stuff must have some impact...