HGridOptimiser Class Reference

#include <HGridOptimiser.h>

Public Member Functions
void	initialise (const MercuryBase &problem, unsigned int numberOfCells, int verbosity)
void	initialisePolyFunc (double omega, std::vector< double > &coeff, unsigned int numberOfCells, int verbosity)
unsigned int	radius2Cell (double r)
	Assigns a BaseParticle of given radius to a certain cell. More...
unsigned int	radius2IntCell (double r)
double	intCell2Min (unsigned int i)
double	intCell2Max (unsigned int i)
double	cell2Min (unsigned int i)
	Computes the left bound of the cell with given ordinal number. More...
double	cell2Max (unsigned int i)
	Computes the right bound of the cell with given ordinal number. More...
double	pdfIntCell (double start, double end, unsigned int i, int p)
double	pdfInt (double start, double end, int power)
	This function calculates: int(f(r)*r^power*dr,r=start..end)/int(f(r)*dr,r=0..omega) with r=a*r+b. More...
double	diffPdfInt (double x, int power)
	diff(int(f(r)*r^power*dr,r=s..e)/int(f(r)*dr,r=0..omega),e)=f(e)*e^power/int(f(r)*dr,r=0..omega) More...
double	expectedCellsIntegralCellNumerator (double start, double end, unsigned int i, int p, double h)
double	diffHExpectedCellsIntegralCellNumerator (double start, double end, unsigned int i, int p, double h)
double	expectedCellsIntegralCellDenominator (double start, double end, unsigned int i)
double	expectedCellsIntegral (double start, double end, int p, double h)
	This function calculates: int((2*r/h+2)^d f(r) dr,r=s..e)/int(f(r) dr,r=s..e)+ Used to calculated the expected number of cells to check at the level with maximum size h for particle radius between start and end. More...
double	diffStartExpectedCellsIntegral (double start, double end, int p, double h)
double	diffEndExpectedCellsIntegral (double start, double end, int p, double h)
double	diffHExpectedCellsIntegral (double start, double end, int p, double h)
double	calculateWork (std::vector< double > &hGridCellSizes, HGridMethod method, int verbosity)
	The amount of work that has to be done to run a simulation using the HGrid, in steps. More...
void	calculateDiffWork (std::vector< double > &hGridCellSizes, std::vector< double > &dfdx, HGridMethod method, int verbosity)
void	calcDfDx (std::vector< double > &hGridCellSizes, std::vector< double > &dfdx, HGridMethod method, int verbosity)
double	checkLimit (std::vector< double > &hGridCellSizes, std::vector< double > &dfdx, int verbosity)
void	applyStep (std::vector< double > &hGridCellSizes, std::vector< double > &dfdx, double stepsize, int verbosity)
double	goldenSectionSearch (std::vector< double > &startHGridCellSizes, std::vector< double > &searchDirection, double min, double cur, double max, HGridMethod method, int verbosity)
void	getOptimalDistribution (std::vector< double > &hGridCellSizes, unsigned int numberOfLevels, HGridMethod method, int verbosity)
void	histNumberParticlesPerCell (std::vector< double > &hGridCellSizes)

Private Attributes
unsigned int	numCells_
	Number of cells, usually called levels in the HGrid. More...
double	rMin_
	Radius of the smallest particle, "rounded" to the largest double that is smaller than the radius of each particle. More...
double	rMax_
	Radius of the largest particle, "rounded" to the smallest double that is larger than the radius of each particle. More...
double	length_
	The weighted length of the domain. More...
double	cellCheckOverContactCheckRatio_
	The ratio of the time required for a single geometric contact detection over the time required to retrieve information from a cell. This seems to be only used for checking the effort required by the HGrid, not to compute the cell sizes. More...
unsigned int	dimension_
	The dimension of the system, usually 3, sometimes 2 or 1. More...
std::vector< double >	cellN_
std::vector< double >	intCellN

Detailed Description

Todo:

Find out what this class does and document it.

Definition at line 36 of file HGridOptimiser.h.

Member Function Documentation

void HGridOptimiser::applyStep	(	std::vector< double > &	hGridCellSizes,
		std::vector< double > &	dfdx,
		double	stepsize,
		int	verbosity
	)

Definition at line 924 of file HGridOptimiser.cc.

Referenced by getOptimalDistribution(), and goldenSectionSearch().

{
    if(verbosity > 0)
    {
        std::cout << "Apply step:" << std::endl;
    }
    for(unsigned int i = 0; i < hGridCellSizes.size() - 1; i++)
    {
        hGridCellSizes[i] += stepsize * dfdx[i];
        if(verbosity > 0)
        {
            std::cout << "hGridCellSizes[i]" << "+" << stepsize << "*" << dfdx[i] << "=" << hGridCellSizes[i] << std::endl;
        }
    }
}

void HGridOptimiser::calcDfDx	(	std::vector< double > &	hGridCellSizes,
		std::vector< double > &	dfdx,
		HGridMethod	method,
		int	verbosity
	)

Definition at line 838 of file HGridOptimiser.cc.

References calculateWork().

{
    if(verbosity > 0)
    {
        std::cout << "calcDfDx" << std::endl;
    }
    double W = calculateWork(hGridCellSizes, method, verbosity - 1);
    double dx = 1e-7;
    double nW;
    dfdx.clear();
    for(unsigned int i = 0; i < hGridCellSizes.size(); i++)
    {
        hGridCellSizes[i] += dx;
        nW = calculateWork(hGridCellSizes, method, verbosity - 1);
        dfdx.push_back((nW - W) / dx);
        hGridCellSizes[i] -= dx;
        if(verbosity > 0)
        {
            std::cout << "i=" << i << " Value=" << hGridCellSizes[i] << " Change=" << nW - W << " dfdx=" << dfdx.back() << std::endl;
        }
    }
}

void HGridOptimiser::calculateDiffWork	(	std::vector< double > &	hGridCellSizes,
		std::vector< double > &	dfdx,
		HGridMethod	method,
		int	verbosity
	)

Definition at line 437 of file HGridOptimiser.cc.

References BOTTOMUP, cellCheckOverContactCheckRatio_, diffEndExpectedCellsIntegral(), diffHExpectedCellsIntegral(), diffPdfInt(), diffStartExpectedCellsIntegral(), dimension_, expectedCellsIntegral(), length_, pdfInt(), and TOPDOWN.

Referenced by getOptimalDistribution().

{
    unsigned int NLevels = hGridCellSizes.size() - 1;
    unsigned int NParams = hGridCellSizes.size();

    if(verbosity > 0)
    {
        std::cout << "Length scale=" << length_ << std::endl;
    }
    if(verbosity > 0)
    {
        std::cout << "Level sizes:";
        for(unsigned int i = 0; i < NParams; i++)
        {
            std::cout << " " << hGridCellSizes[i];
        }
        std::cout << std::endl;
    }

    std::vector<double> particlesPerLevel;
    for(unsigned int i = 0; i < NLevels; i++)
    {
        particlesPerLevel.push_back(pdfInt(0.5 * hGridCellSizes[i], 0.5 * hGridCellSizes[i + 1], 0));
        if(particlesPerLevel.back() < 0)
        {
            particlesPerLevel.back() = 0;
        }
    }
    if(verbosity > 0)
    {
        std::cout << "Particles per level: ";
        for(unsigned int i = 0; i < NLevels; i++)
        {
            std::cout << " " << particlesPerLevel[i];
        }
        std::cout << std::endl;
    }

    //How changes particlesPerLeve[i] when hGridCellSize[j] is changed
    std::vector<std::vector<double> > diffParticlesPerLevel;
    diffParticlesPerLevel.resize(NLevels);
    for(unsigned int i = 0; i < NLevels; i++)
    {
        for(unsigned int j = 0; j < NParams; j++)
        {
            diffParticlesPerLevel[i].push_back(0.0);
            if(j == i)
            {
                diffParticlesPerLevel[i][j] += -0.5 * diffPdfInt(0.5 * hGridCellSizes[i], 0);
            }
            if(j == i + 1)
            {
                diffParticlesPerLevel[i][j] += 0.5 * diffPdfInt(0.5 * hGridCellSizes[i + 1], 0);
            }
        }
    }
    if(verbosity > 0)
    {
        std::cout << "Diff Particles per level: " << std::endl;
        for(unsigned int i = 0; i < NLevels; i++)
        {
            for(unsigned int j = 0; j < NParams; j++)
            {
                std::cout << " " << std::setw(12) << diffParticlesPerLevel[i][j];
            }
            std::cout << std::endl;
        }
    }

    std::vector<double> cellsPerLevel;
    for(unsigned int i = 0; i < NLevels; i++)
    {
        cellsPerLevel.push_back(pow(length_ / hGridCellSizes[i + 1], dimension_));
    }
    if(verbosity > 0)
    {
        std::cout << "Cells per level: ";
        for(unsigned int i = 0; i < NLevels; i++)
        {
            std::cout << " " << cellsPerLevel[i];
        }
        std::cout << std::endl;
    }

    //How changes cellsPerLevel[i] when hGridCellSize[j] is changed    
    std::vector<std::vector<double> > diffCellsPerLevel;
    diffCellsPerLevel.resize(hGridCellSizes.size());
    for(unsigned int i = 0; i < NLevels; i++)
    {
        for(unsigned int j = 0; j < NParams; j++)
        {
            if(j == i + 1)
            {
                diffCellsPerLevel[i].push_back(-pow(length_ / hGridCellSizes[i + 1], dimension_) * dimension_ / hGridCellSizes[i + 1]);
            }
            else
            {
                diffCellsPerLevel[i].push_back(0.0);
            }
        }
    }
    if(verbosity > 0)
    {
        std::cout << "Diff Cells per level: " << std::endl;
        for(unsigned int i = 0; i < NLevels; i++)
        {
            for(unsigned int j = 0; j < NParams; j++)
            {
                std::cout << " " << std::setw(12) << diffCellsPerLevel[i][j];
            }
            std::cout << std::endl;
        }
    }

    std::vector<double> particlesPerCell;
    for(unsigned int i = 0; i < NLevels; i++)
    {
        particlesPerCell.push_back(particlesPerLevel[i] / cellsPerLevel[i]);
    }
    if(verbosity > 0)
    {
        std::cout << "Particles per cell: ";
        for(unsigned int i = 0; i < particlesPerCell.size(); i++)
        {
            std::cout << " " << particlesPerCell[i];
        }
        std::cout << std::endl;
    }

    //How changes particlesPerCell[i] when hGridCellSize[j] is changed    
    std::vector<std::vector<double> > diffParticlesPerCell;
    diffParticlesPerCell.resize(NLevels);
    for(unsigned int i = 0; i < NLevels; i++)
    {
        for(unsigned int j = 0; j < NParams; j++)
        {
            diffParticlesPerCell[i].push_back((diffParticlesPerLevel[i][j] * cellsPerLevel[i] - particlesPerLevel[i] * diffCellsPerLevel[i][j]) / (cellsPerLevel[i] * cellsPerLevel[i]));
        }
    }
    if(verbosity > 0)
    {
        std::cout << "Diff Particles per Cell: " << std::endl;
        for(unsigned int i = 0; i < NLevels; i++)
        {
            for(unsigned int j = 0; j < NParams; j++)
            {
                std::cout << " " << std::setw(12) << diffParticlesPerCell[i][j];
            }
            std::cout << std::endl;
        }
    }

    double contactWork = 0, overheadWork = 0;
    std::vector<double> diffContactWork;
    std::vector<double> diffOverheadWork;
    diffContactWork.resize(NParams);
    diffOverheadWork.resize(NParams);
    for(unsigned int j = 0; j < NParams; j++)
    {
        diffContactWork[j] = 0;
        diffOverheadWork[j] = 0;
    }

    for(unsigned int i = 0; i < NLevels; i++)
    {
        contactWork += 0.5 * pow(3, dimension_) * particlesPerCell[i] * particlesPerLevel[i];
        overheadWork += 0.5 * (pow(3, dimension_) + 1) * particlesPerLevel[i];
        for(unsigned int j = 0; j < NParams; j++)
        {
            diffContactWork[j] += 0.5 * pow(3, dimension_) * (diffParticlesPerCell[i][j] * particlesPerLevel[i] + particlesPerCell[i] * diffParticlesPerLevel[i][j]);
            diffOverheadWork[j] += 0.5 * (pow(3, dimension_) + 1) * diffParticlesPerLevel[i][j];
        }

        unsigned int startJ, endJ;
        switch(method)
        {
        case BOTTOMUP:
        {
            startJ = i + 1;
            endJ = NLevels;
            break;
        }
        case TOPDOWN:
        {
            startJ = 0;
            endJ = i;
            break;
        }
        }
        for(unsigned int j = startJ; j < endJ; j++)
        {
            double numberOfCellToCheck;
            numberOfCellToCheck = expectedCellsIntegral(0.5 * hGridCellSizes[i], 0.5 * hGridCellSizes[i + 1], dimension_, hGridCellSizes[j + 1]);

            std::vector<double> diffNumberOfCellToCheck;
            for(unsigned int k = 0; k < NParams; k++)
            {
                diffNumberOfCellToCheck.push_back(0.0);
                if(k == i)
                {
                    diffNumberOfCellToCheck[k] += 0.5 * diffStartExpectedCellsIntegral(0.5 * hGridCellSizes[i], 0.5 * hGridCellSizes[i + 1], dimension_, hGridCellSizes[j + 1]);
                }
                if(k == i + 1)
                {
                    diffNumberOfCellToCheck[k] += 0.5 * diffEndExpectedCellsIntegral(0.5 * hGridCellSizes[i], 0.5 * hGridCellSizes[i + 1], dimension_, hGridCellSizes[j + 1]);
                }
                if(k == j + 1)
                {
                    diffNumberOfCellToCheck[k] += diffHExpectedCellsIntegral(0.5 * hGridCellSizes[i], 0.5 * hGridCellSizes[i + 1], dimension_, hGridCellSizes[j + 1]);
                }
            }
            contactWork += particlesPerLevel[i] * numberOfCellToCheck * particlesPerCell[j];
            overheadWork += particlesPerLevel[i] * numberOfCellToCheck;
            //std::cout<<"i= "<<i<<" j= "<<j<<" NumberOfCellToCheck= "<<numberOfCellToCheck<<" diffNumberOfCellToCheck[2]= "<<diffNumberOfCellToCheck[2]<<std::endl;
            for(unsigned int k = 0; k < NParams; k++)
            {
                diffContactWork[k] += (diffParticlesPerLevel[i][k] * numberOfCellToCheck * particlesPerCell[j] + particlesPerLevel[i] * diffNumberOfCellToCheck[k] * particlesPerCell[j] + particlesPerLevel[i] * numberOfCellToCheck * diffParticlesPerCell[j][k]);
                diffOverheadWork[k] += (diffParticlesPerLevel[i][k] * numberOfCellToCheck + particlesPerLevel[i] * diffNumberOfCellToCheck[k]);
            }
        }
    }
    if(verbosity > 0)
    {
        std::cout << "Contact work: " << contactWork << std::endl;
        std::cout << "Overhead work: " << cellCheckOverContactCheckRatio_ * overheadWork << std::endl;
        std::cout << "Sum work: " << contactWork + cellCheckOverContactCheckRatio_ * overheadWork << std::endl;
        std::cout << "diff Contactwork: ";
        for(unsigned int j = 0; j < NParams; j++)
        {
            std::cout << " " << diffContactWork[j];
        }
        std::cout << std::endl;
        std::cout << "diff Overheadwork: ";
        for(unsigned int j = 0; j < NParams; j++)
        {
            std::cout << " " << cellCheckOverContactCheckRatio_ * diffOverheadWork[j];
        }
        std::cout << std::endl;
        std::cout << "diff sum: ";
        for(unsigned int j = 0; j < NParams; j++)
        {
            std::cout << " " << diffContactWork[j] + cellCheckOverContactCheckRatio_ * diffOverheadWork[j];
        }
        std::cout << std::endl;
    }
    dfdx.clear();
    for(unsigned int j = 0; j < NParams; j++)
    {
        dfdx.push_back(diffContactWork[j] + cellCheckOverContactCheckRatio_ * diffOverheadWork[j]);
    }
}

double HGridOptimiser::calculateWork	(	std::vector< double > &	hGridCellSizes,
		HGridMethod	method,
		int	verbosity
	)

The amount of work that has to be done to run a simulation using the HGrid, in steps.

Definition at line 689 of file HGridOptimiser.cc.

References BOTTOMUP, cellCheckOverContactCheckRatio_, dimension_, expectedCellsIntegral(), length_, pdfInt(), and TOPDOWN.

Referenced by calcDfDx(), getOptimalDistribution(), and goldenSectionSearch().

{
    unsigned int NLevels = hGridCellSizes.size() - 1;
    unsigned int NParams = hGridCellSizes.size();

    if(verbosity > 1)
    {
        std::cout << "Length scale=" << length_ << std::endl;
    }
    if(verbosity > 0)
    {
        std::cout << "Level sizes:" << std::endl;
        for(unsigned int i = 0; i < NParams; i++)
        {
            std::cout << std::setw(10) << hGridCellSizes[i] << " ";
        }
        std::cout << std::endl;
    }

    std::vector<double> particlesPerLevel;
    for(unsigned int i = 0; i < NLevels; i++)
    {
        particlesPerLevel.push_back(pdfInt(0.5 * hGridCellSizes[i], 0.5 * hGridCellSizes[i + 1], 0));
        if(particlesPerLevel.back() < 0)
        {
            particlesPerLevel.back() = 0;
        }
    }
    if(verbosity > 0)
    {
        std::cout << "Particles per level:" << std::endl;
        for(unsigned int i = 0; i < NLevels; i++)
        {
            std::cout << std::setw(10) << particlesPerLevel[i] << " ";
        }
        std::cout << std::endl;
    }

    std::vector<double> cellsPerLevel;
    for(unsigned int i = 0; i < NLevels; i++)
    {
        cellsPerLevel.push_back(pow(length_ / hGridCellSizes[i + 1], dimension_));
    }
    if(verbosity > 1)
    {
        std::cout << "Cells per level:" << std::endl;
        for(unsigned int i = 0; i < NLevels; i++)
        {
            std::cout << std::setw(10) << cellsPerLevel[i] << " ";
        }
        std::cout << std::endl;
    }

    std::vector<double> particlesPerCell;
    for(unsigned int i = 0; i < NLevels; i++)
    {
        particlesPerCell.push_back(particlesPerLevel[i] / cellsPerLevel[i]);
    }
    if(verbosity > 1)
    {
        std::cout << "Particles per cell:" << std::endl;
        for(unsigned int i = 0; i < particlesPerCell.size(); i++)
        {
            std::cout << std::setw(10) << particlesPerCell[i] << " ";
        }
        std::cout << std::endl;
    }

    std::vector<std::vector<double> > contactWork;
    std::vector<std::vector<double> > overheadWork;
    contactWork.resize(NLevels);
    overheadWork.resize(NLevels);
    for(unsigned int i = 0; i < NLevels; i++)
    {
        contactWork[i].resize(NLevels);
        overheadWork[i].resize(NLevels);
        for(unsigned int j = 0; j < NLevels; j++)
        {
            contactWork[i][j] = 0;
            overheadWork[i][j] = 0;
        }
    }

    for(unsigned int i = 0; i < NLevels; i++)
    {
        contactWork[i][i] += 0.5 * pow(3, dimension_) * particlesPerCell[i] * particlesPerLevel[i];
        overheadWork[i][i] += 0.5 * (pow(3, dimension_) + 1) * particlesPerLevel[i];

        unsigned int startJ, endJ;
        switch(method)
        {
        case BOTTOMUP:
        {
            startJ = i + 1;
            endJ = NLevels;
            break;
        }
        case TOPDOWN:
        {
            startJ = 0;
            endJ = i;
            break;
        }
        }
        for(unsigned int j = startJ; j < endJ; j++)
        {
            double numberOfCellToCheck;
            numberOfCellToCheck = expectedCellsIntegral(0.5 * hGridCellSizes[i], 0.5 * hGridCellSizes[i + 1], dimension_, hGridCellSizes[j + 1]);
            contactWork[i][j] += particlesPerLevel[i] * numberOfCellToCheck * particlesPerCell[j];
            overheadWork[i][j] += particlesPerLevel[i] * numberOfCellToCheck;

        }
    }

    double totalContactWork = 0, totalOverheadWork = 0;
    for(unsigned int i = 0; i < NLevels; i++)
    {
        for(unsigned int j = 0; j < NLevels; j++)
        {
            totalContactWork += contactWork[i][j];
            totalOverheadWork += overheadWork[i][j];
        }
    }

    if(verbosity > 0)
    {
        std::cout << "Contact work: " << totalContactWork << std::endl;
        for(unsigned int i = 0; i < NLevels; i++)
        {
            for(unsigned int j = 0; j < NLevels; j++)
            {
                std::cout << std::setw(10) << contactWork[i][j] << " ";
            }
            std::cout << std::endl;
        }
        std::cout << "Overhead work: " << totalOverheadWork << std::endl;
        for(unsigned int i = 0; i < NLevels; i++)
        {
            for(unsigned int j = 0; j < NLevels; j++)
            {
                std::cout << std::setw(10) << overheadWork[i][j] << " ";
            }
            std::cout << std::endl;
        }
        std::cout << "Total work: " << totalContactWork + totalOverheadWork << std::endl;
    }
    return totalContactWork + cellCheckOverContactCheckRatio_ * totalOverheadWork;
}

double HGridOptimiser::cell2Max

(

unsigned int

i

)

Computes the right bound of the cell with given ordinal number.

Computes the right bound of the cell with given ordinal number.

Definition at line 199 of file HGridOptimiser.cc.

References numCells_, rMax_, and rMin_.

Referenced by initialise(), and initialisePolyFunc().

{
    return rMin_ + (rMax_ - rMin_) * (i + 1) / numCells_;
}

double HGridOptimiser::cell2Min

(

unsigned int

i

)

Computes the left bound of the cell with given ordinal number.

Computes the left bound of the cell with given ordinal number.

Definition at line 191 of file HGridOptimiser.cc.

References numCells_, rMax_, and rMin_.

Referenced by initialise(), and initialisePolyFunc().

{
    return rMin_ + (rMax_ - rMin_) * i / numCells_;
}

double HGridOptimiser::checkLimit	(	std::vector< double > &	hGridCellSizes,
		std::vector< double > &	dfdx,
		int	verbosity
	)

Definition at line 861 of file HGridOptimiser.cc.

Referenced by getOptimalDistribution().

{
    if(verbosity > 0)
    {
        std::cout << "checkLimits part 1 remove levels" << std::endl;
    }

    dfdx[0] = 0;
    dfdx.back() = 0;
    for(unsigned int i = 1; i < hGridCellSizes.size(); i++)
    {
        if(std::abs(hGridCellSizes[i - 1] - hGridCellSizes[i]) < 1e-6)
        {
            if(dfdx[i - 1] < dfdx[i])
            {
                if(verbosity > 0)
                {
                    std::cout << "Remove level" << i << std::endl;
                }
                if(i > 0.5 * hGridCellSizes.size())
                {
                    hGridCellSizes.erase(hGridCellSizes.begin() + i - 1);
                    dfdx.erase(dfdx.begin() + i - 1);
                }
                else
                {
                    hGridCellSizes.erase(hGridCellSizes.begin() + i);
                    dfdx.erase(dfdx.begin() + i);
                }
                i--;
            }
        }
    }

    if(verbosity > 0)
    {
        std::cout << "checkLimits part 2 calculate limit" << std::endl;
        for(unsigned int i = 0; i < hGridCellSizes.size(); i++)
        {
            std::cout << "i=" << i << " Value=" << hGridCellSizes[i] << " dfdx=" << dfdx[i] << std::endl;
        }
    }
    double maxStepSize = -std::numeric_limits<double>::max();
    double nmax;
    for(unsigned int i = 1; i < hGridCellSizes.size(); i++)
    {
        nmax = (hGridCellSizes[i - 1] - hGridCellSizes[i]) / (dfdx[i] - dfdx[i - 1]);
        if(verbosity > 0)
        {
            std::cout << "Max f=" << nmax << " because otherwise " << hGridCellSizes[i - 1] << "+x*" << dfdx[i - 1] << ">" << hGridCellSizes[i] << "+x*" << dfdx[i] << std::endl;
        }
        if(nmax < 0)
        {
            maxStepSize = std::max(nmax, maxStepSize);
        }
    }
    if(verbosity > 0)
    {
        std::cout << "maxStepSize=" << maxStepSize << std::endl;
    }
    return maxStepSize;
}

double HGridOptimiser::diffEndExpectedCellsIntegral	(	double	start,
		double	end,
		int	p,
		double	h
	)

Definition at line 372 of file HGridOptimiser.cc.

References expectedCellsIntegralCellDenominator(), expectedCellsIntegralCellNumerator(), intCell2Max(), intCell2Min(), intCellN, and radius2IntCell().

Referenced by calculateDiffWork().

{
    unsigned int startCell = radius2IntCell(start);
    unsigned int endCell = radius2IntCell(end);

    double numerator = 0;
    double denominator = 0;

    double a = (intCellN[endCell + 1] - intCellN[endCell]) / (intCell2Max(endCell) - intCell2Min(endCell));
    double b = (intCellN[endCell] * intCell2Max(endCell) - intCellN[endCell + 1] * intCell2Min(endCell)) / (intCell2Max(endCell) - intCell2Min(endCell));

    double diffTeller = (2.0 / h * p * std::pow(2.0 * (end + h) / h, p - 1) * (a * p * end - a * h + a * end + b * p + 2 * b) * (end + h) / ((p + 1) * (p + 2)) +
                         std::pow(2.0 * (end + h) / h, p) * (a * p + a) * (end + h) / ((p + 1) * (p + 2)) +
                         std::pow(2.0 * (end + h) / h, p) * (a * p * end - a * h + a * end + b * p + 2 * b) / ((p + 1) * (p + 2)));
    double diffNoemer = (a * end + b);
    if(startCell == endCell)
    {
        numerator = expectedCellsIntegralCellNumerator(start, end, startCell, p, h);
        denominator = expectedCellsIntegralCellDenominator(start, end, startCell);
    }
    else
    {
        numerator = expectedCellsIntegralCellNumerator(start, intCell2Max(startCell), startCell, p, h);
        denominator = expectedCellsIntegralCellDenominator(start, intCell2Max(startCell), startCell);
        for(unsigned int i = startCell + 1; i < endCell; i++)
        {
            numerator += expectedCellsIntegralCellNumerator(intCell2Min(i), intCell2Max(i), i, p, h);
            denominator += expectedCellsIntegralCellDenominator(intCell2Min(i), intCell2Max(i), i);
        }
        numerator += expectedCellsIntegralCellNumerator(intCell2Min(endCell), end, endCell, p, h);
        denominator += expectedCellsIntegralCellDenominator(intCell2Min(endCell), end, endCell);

    }
    //std::cout<<"new: numerator="<<numerator<<" denominator="<<denominator<<" diffTeller="<<diffTeller<<" diffNoemer="<<diffNoemer<<std::endl;
    return(diffTeller * denominator - numerator * diffNoemer) / (denominator * denominator);
}

double HGridOptimiser::diffHExpectedCellsIntegral	(	double	start,
		double	end,
		int	p,
		double	h
	)

Definition at line 409 of file HGridOptimiser.cc.

References diffHExpectedCellsIntegralCellNumerator(), expectedCellsIntegralCellDenominator(), intCell2Max(), intCell2Min(), and radius2IntCell().

Referenced by calculateDiffWork().

{
    unsigned int startCell = radius2IntCell(start);
    unsigned int endCell = radius2IntCell(end);

    double denominator = 0;
    double diffTeller = 0;

    if(startCell == endCell)
    {
        diffTeller = diffHExpectedCellsIntegralCellNumerator(start, end, startCell, p, h);
        denominator = expectedCellsIntegralCellDenominator(start, end, startCell);
    }
    else
    {
        diffTeller = diffHExpectedCellsIntegralCellNumerator(start, intCell2Max(startCell), startCell, p, h);
        denominator = expectedCellsIntegralCellDenominator(start, intCell2Max(startCell), startCell);
        for(unsigned int i = startCell + 1; i < endCell; i++)
        {
            diffTeller += diffHExpectedCellsIntegralCellNumerator(intCell2Min(i), intCell2Max(i), i, p, h);
            denominator += expectedCellsIntegralCellDenominator(intCell2Min(i), intCell2Max(i), i);
        }
        diffTeller += diffHExpectedCellsIntegralCellNumerator(intCell2Min(endCell), end, endCell, p, h);
        denominator += expectedCellsIntegralCellDenominator(intCell2Min(endCell), end, endCell);
    }
    return diffTeller / denominator;
}

double HGridOptimiser::diffHExpectedCellsIntegralCellNumerator	(	double	start,
		double	end,
		unsigned int	i,
		int	p,
		double	h
	)

Definition at line 276 of file HGridOptimiser.cc.

References intCell2Max(), intCell2Min(), and intCellN.

Referenced by diffHExpectedCellsIntegral().

{
    double a = (intCellN[i + 1] - intCellN[i]) / (intCell2Max(i) - intCell2Min(i));
    double b = (intCellN[i] * intCell2Max(i) - intCellN[i + 1] * intCell2Min(i)) / (intCell2Max(i) - intCell2Min(i));
    double r;
    r = start;
    //double min =std::pow(2.0*r/h+2.0,p)*(a*p*r-a*h+a*r+b*p+2*b)*(r+h)/((p+1)*(p+2));
    double min = (-2.0 * r / h / h * p * std::pow(2.0 * r / h + 2.0, p - 1) * (a * p * r - a * h + a * r + b * p + 2 * b) * (r + h) +
                  -std::pow(2.0 * r / h + 2.0, p) * a * (r + h) +
                  std::pow(2.0 * r / h + 2.0, p) * (a * p * r - a * h + a * r + b * p + 2 * b)) / ((p + 1) * (p + 2));
    r = end;
    //double plus=std::pow(2.0*r/h+2.0,p)*(a*p*r-a*h+a*r+b*p+2*b)*(r+h)/((p+1)*(p+2));
    double plus = (-2.0 * r / h / h * p * std::pow(2.0 * r / h + 2.0, p - 1) * (a * p * r - a * h + a * r + b * p + 2 * b) * (r + h) +
                   -std::pow(2.0 * r / h + 2.0, p) * a * (r + h) +
                   std::pow(2.0 * r / h + 2.0, p) * (a * p * r - a * h + a * r + b * p + 2 * b)) / ((p + 1) * (p + 2));
    return plus - min;
}

double HGridOptimiser::diffPdfInt	(	double	x,
		int	power
	)

diff(int(f(r)*r^power*dr,r=s..e)/int(f(r)*dr,r=0..omega),e)=f(e)*e^power/int(f(r)*dr,r=0..omega)

Definition at line 251 of file HGridOptimiser.cc.

References intCell2Max(), intCell2Min(), intCellN, numCells_, pdfIntCell(), and radius2IntCell().

Referenced by calculateDiffWork().

{
    unsigned int xCell = radius2IntCell(x);
    double denominator = 0;
    for(unsigned int i = 0; i <= numCells_; i++)
    {
        denominator += pdfIntCell(intCell2Min(i), intCell2Max(i), i, 0);
    }
    double a = (intCellN[xCell + 1] - intCellN[xCell]) / (intCell2Max(xCell) - intCell2Min(xCell));
    double b = (intCellN[xCell] * intCell2Max(xCell) - intCellN[xCell + 1] * intCell2Min(xCell)) / (intCell2Max(xCell) - intCell2Min(xCell));
    return(a * std::pow(x, p + 1) + b * std::pow(x, p)) / denominator;
}

double HGridOptimiser::diffStartExpectedCellsIntegral	(	double	start,
		double	end,
		int	p,
		double	h
	)

Definition at line 335 of file HGridOptimiser.cc.

References expectedCellsIntegralCellDenominator(), expectedCellsIntegralCellNumerator(), intCell2Max(), intCell2Min(), intCellN, and radius2IntCell().

Referenced by calculateDiffWork().

{
    unsigned int startCell = radius2IntCell(start);
    unsigned int endCell = radius2IntCell(end);

    double numerator = 0;
    double denominator = 0;

    double a = (intCellN[startCell + 1] - intCellN[startCell]) / (intCell2Max(startCell) - intCell2Min(startCell));
    double b = (intCellN[startCell] * intCell2Max(startCell) - intCellN[startCell + 1] * intCell2Min(startCell)) / (intCell2Max(startCell) - intCell2Min(startCell));

    double diffTeller = -(2.0 / h * p * std::pow(2.0 * (start + h) / h, p - 1) * (a * p * start - a * h + a * start + b * p + 2 * b) * (start + h) / ((p + 1) * (p + 2)) +
                          std::pow(2.0 * (start + h) / h, p) * (a * p + a) * (start + h) / ((p + 1) * (p + 2)) +
                          std::pow(2.0 * (start + h) / h, p) * (a * p * start - a * h + a * start + b * p + 2 * b) / ((p + 1) * (p + 2)));
    double diffNoemer = -(a * start + b);
    if(startCell == endCell)
    {
        numerator = expectedCellsIntegralCellNumerator(start, end, startCell, p, h);
        denominator = expectedCellsIntegralCellDenominator(start, end, startCell);
    }
    else
    {
        numerator = expectedCellsIntegralCellNumerator(start, intCell2Max(startCell), startCell, p, h);
        denominator = expectedCellsIntegralCellDenominator(start, intCell2Max(startCell), startCell);
        for(unsigned int i = startCell + 1; i < endCell; i++)
        {
            numerator += expectedCellsIntegralCellNumerator(intCell2Min(i), intCell2Max(i), i, p, h);
            denominator += expectedCellsIntegralCellDenominator(intCell2Min(i), intCell2Max(i), i);
        }
        numerator += expectedCellsIntegralCellNumerator(intCell2Min(endCell), end, endCell, p, h);
        denominator += expectedCellsIntegralCellDenominator(intCell2Min(endCell), end, endCell);

    }
    //std::cout<<"new: numerator="<<numerator<<" denominator="<<denominator<<" diffTeller="<<diffTeller<<" diffNoemer="<<diffNoemer<<std::endl;
    return(diffTeller * denominator - numerator * diffNoemer) / (denominator * denominator);
}

double HGridOptimiser::expectedCellsIntegral	(	double	start,
		double	end,
		int	p,
		double	h
	)

This function calculates: int((2*r/h+2)^d f(r) dr,r=s..e)/int(f(r) dr,r=s..e)+ Used to calculated the expected number of cells to check at the level with maximum size h for particle radius between start and end.

Definition at line 305 of file HGridOptimiser.cc.

References expectedCellsIntegralCellDenominator(), expectedCellsIntegralCellNumerator(), intCell2Max(), intCell2Min(), and radius2IntCell().

Referenced by calculateDiffWork(), and calculateWork().

{
    unsigned int startCell = radius2IntCell(start);
    unsigned int endCell = radius2IntCell(end);

    double numerator = 0;
    double denominator = 0;
    if(startCell == endCell)
    {
        numerator = expectedCellsIntegralCellNumerator(start, end, startCell, p, h);
        denominator = expectedCellsIntegralCellDenominator(start, end, startCell);
    }
    else
    {
        numerator = expectedCellsIntegralCellNumerator(start, intCell2Max(startCell), startCell, p, h);
        denominator = expectedCellsIntegralCellDenominator(start, intCell2Max(startCell), startCell);
        for(unsigned int i = startCell + 1; i < endCell; i++)
        {
            numerator += expectedCellsIntegralCellNumerator(intCell2Min(i), intCell2Max(i), i, p, h);
            denominator += expectedCellsIntegralCellDenominator(intCell2Min(i), intCell2Max(i), i);
        }
        numerator += expectedCellsIntegralCellNumerator(intCell2Min(endCell), end, endCell, p, h);
        denominator += expectedCellsIntegralCellDenominator(intCell2Min(endCell), end, endCell);

    }
    //    std::cout<<"New: numerator="<<numerator<<" denominator="<<denominator<<std::endl;
    return numerator / denominator;
}

double HGridOptimiser::expectedCellsIntegralCellDenominator	(	double	start,
		double	end,
		unsigned int	i
	)

Definition at line 294 of file HGridOptimiser.cc.

References intCell2Max(), intCell2Min(), and intCellN.

Referenced by diffEndExpectedCellsIntegral(), diffHExpectedCellsIntegral(), diffStartExpectedCellsIntegral(), and expectedCellsIntegral().

{
    double a = (intCellN[i + 1] - intCellN[i]) / (intCell2Max(i) - intCell2Min(i));
    double b = (intCellN[i] * intCell2Max(i) - intCellN[i + 1] * intCell2Min(i)) / (intCell2Max(i) - intCell2Min(i));
    return(a / 2.0 * (pow(end, 2) - pow(start, 2)) + b * (pow(end, 1) - pow(start, 1)));
}

double HGridOptimiser::expectedCellsIntegralCellNumerator	(	double	start,
		double	end,
		unsigned int	i,
		int	p,
		double	h
	)

Definition at line 264 of file HGridOptimiser.cc.

References intCell2Max(), intCell2Min(), and intCellN.

Referenced by diffEndExpectedCellsIntegral(), diffStartExpectedCellsIntegral(), and expectedCellsIntegral().

{
    double a = (intCellN[i + 1] - intCellN[i]) / (intCell2Max(i) - intCell2Min(i));
    double b = (intCellN[i] * intCell2Max(i) - intCellN[i + 1] * intCell2Min(i)) / (intCell2Max(i) - intCell2Min(i));
    double r;
    r = start;
    double min = std::pow(2.0 * (r + h) / h, p) * (a * p * r - a * h + a * r + b * p + 2 * b) * (r + h) / ((p + 1) * (p + 2));
    r = end;
    double plus = std::pow(2.0 * (r + h) / h, p) * (a * p * r - a * h + a * r + b * p + 2 * b) * (r + h) / ((p + 1) * (p + 2));
    return plus - min;
}

void HGridOptimiser::getOptimalDistribution	(	std::vector< double > &	hGridCellSizes,
		unsigned int	numberOfLevels,
		HGridMethod	method,
		int	verbosity
	)

Definition at line 998 of file HGridOptimiser.cc.

References applyStep(), calculateDiffWork(), calculateWork(), checkLimit(), goldenSectionSearch(), rMax_, and rMin_.

{
    hGridCellSizes.clear();
    for(unsigned int i = 0; i <= numberOfLevels; i++)
    {
        hGridCellSizes.push_back(2.0 * (rMin_ + (rMax_ - rMin_) * i / numberOfLevels));
    }

    std::vector<double> dfdx;
    double step, max, W, oW;
    W = calculateWork(hGridCellSizes, method, verbosity - 3);
    int stepNumber = 0;
    do
    {
        oW = W;
        calculateDiffWork(hGridCellSizes, dfdx, method, verbosity - 3);
        max = checkLimit(hGridCellSizes, dfdx, verbosity - 3);
        step = goldenSectionSearch(hGridCellSizes, dfdx, max, 0.5 * max, 0, method, verbosity - 2);
        applyStep(hGridCellSizes, dfdx, step, verbosity - 3);
        W = calculateWork(hGridCellSizes, method, verbosity - 3);
        while(W > oW)
        {
            if(verbosity > 1)
            {
                std::cout << stepNumber << " Bad step, trying closer search range" << std::endl;
            }
            applyStep(hGridCellSizes, dfdx, -step, verbosity - 3);
            max /= 2;
            step = goldenSectionSearch(hGridCellSizes, dfdx, max, 0.5 * max, 0, method, verbosity - 2);
            applyStep(hGridCellSizes, dfdx, step, verbosity - 3);
            W = calculateWork(hGridCellSizes, method, verbosity - 3);
        }
        ++stepNumber;
        if(verbosity > 1)
        {
            std::cout << stepNumber << " Correct step: old work=" << oW << " new work=" << W << " difference=" << oW - W << " step=" << step / max << std::endl;
        }
    }
    while(oW - W > 1e-13 && stepNumber);
    calculateDiffWork(hGridCellSizes, dfdx, method, verbosity - 2);
    if(verbosity > 0)
    {
        std::cout << "Work=" << W << std::endl;
        std::cout << "Optimal cell sizes:";
        for(unsigned int i = 0; i < hGridCellSizes.size(); i++)
        {
            std::cout << " " << hGridCellSizes[i];
        }
        std::cout << std::endl;
    }
}

double HGridOptimiser::goldenSectionSearch	(	std::vector< double > &	startHGridCellSizes,
		std::vector< double > &	searchDirection,
		double	min,
		double	cur,
		double	max,
		HGridMethod	method,
		int	verbosity
	)

Definition at line 940 of file HGridOptimiser.cc.

References applyStep(), and calculateWork().

Referenced by getOptimalDistribution().

{
    std::vector<double> curHGridCellSizes = startHGridCellSizes;
    std::vector<double> xHGridCellSizes = startHGridCellSizes;

    applyStep(curHGridCellSizes, searchDirection, cur, verbosity - 1);
    double curWork = calculateWork(curHGridCellSizes, method, verbosity - 1);

    double x;
    double resphi = 2.0 - 0.5 * (1 + std::sqrt(5));

    //Determine new probing point
    if(max - cur > cur - min)
    {
        //x between cur and max
        x = cur + resphi * (max - cur);
    }
    else
    {
        //x between min and cur
        x = cur - resphi * (cur - min);
    }

    if(std::abs(max - min) < 1e-7 * (std::abs(cur) + std::abs(x)))
    {
        return 0.5 * (min + max);
    }

    applyStep(xHGridCellSizes, searchDirection, x, 0);
    double xWork = calculateWork(xHGridCellSizes, method, 0);
    if(verbosity > 0)
    {
        std::cout << "min=" << min << " max=" << max << " cur=" << cur << " curWork=" << curWork << " x=" << x << " xWork=" << xWork << std::endl;
    }
    if(xWork < curWork) //x is better
    {
        if(max - cur > cur - min)
        {
            return goldenSectionSearch(startHGridCellSizes, searchDirection, cur, x, max, method, verbosity);
        }
        else
        {
            return goldenSectionSearch(startHGridCellSizes, searchDirection, min, x, cur, method, verbosity);
        }
    }
    else //cur is better
    {
        if(max - cur > cur - min)
        {
            return goldenSectionSearch(startHGridCellSizes, searchDirection, min, cur, x, method, verbosity);
        }
        else
        {
            return goldenSectionSearch(startHGridCellSizes, searchDirection, x, cur, max, method, verbosity);
        }
    }
}

void HGridOptimiser::histNumberParticlesPerCell

(

std::vector< double > &

hGridCellSizes

)

Definition at line 1050 of file HGridOptimiser.cc.

References dimension_, mathsFunc::factorial(), length_, and pdfInt().

{
    unsigned int NLevels = hGridCellSizes.size() - 1;

    std::vector<double> particlesPerLevel;
    for(unsigned int i = 0; i < NLevels; i++)
    {
        particlesPerLevel.push_back(pdfInt(0.5 * hGridCellSizes[i], 0.5 * hGridCellSizes[i + 1], 0));
        if(particlesPerLevel.back() < 0)
        {
            particlesPerLevel.back() = 0;
        }
    }

    std::vector<double> cellsPerLevel;
    for(unsigned int i = 0; i < NLevels; i++)
    {
        cellsPerLevel.push_back(pow(length_ / hGridCellSizes[i + 1], dimension_));
    }

    for(unsigned int i = 0; i < NLevels; i++)
    {
        double l = particlesPerLevel[i] / cellsPerLevel[i];
        std::cout << "Histogram for level " << i << ":";
        for(unsigned int k = 0; k <= 10; k++)
        {
            std::cout << " " << std::setw(8) << std::floor(std::pow(l, k) * std::exp(-l) / mathsFunc::factorial(k)*1e4 * cellsPerLevel[i]);
        }
        std::cout << std::endl;
    }
}

void HGridOptimiser::initialise	(	const MercuryBase &	problem,
		unsigned int	numberOfCells,
		int	verbosity
	)

Definition at line 33 of file HGridOptimiser.cc.

References BaseHandler< T >::begin(), cell2Max(), cell2Min(), cellCheckOverContactCheckRatio_, cellN_, dimension_, BaseHandler< T >::end(), BaseParticle::getInteractionRadius(), ParticleHandler::getLargestParticle(), BaseHandler< T >::getNumberOfObjects(), ParticleHandler::getSmallestParticle(), DPMBase::getSystemDimensions(), DPMBase::getXMax(), DPMBase::getXMin(), DPMBase::getYMax(), DPMBase::getYMin(), DPMBase::getZMax(), DPMBase::getZMin(), intCell2Max(), intCell2Min(), intCellN, length_, numCells_, DPMBase::particleHandler, radius2Cell(), rMax_, and rMin_.

{
    //assign the number of cells
    numCells_ = numberOfCells;
    cellCheckOverContactCheckRatio_ = 0.2;
    
    cellN_.resize(numCells_);
    for(std::vector<double>::iterator it = cellN_.begin(); it != cellN_.end(); ++it)
    {
        *it = 0.0;
    }

    //Set the minimum and maximum radius of the particles.
    rMin_ = nextafter(problem.particleHandler.getSmallestParticle()->getInteractionRadius(), 0.0);
    rMax_ = nextafter(nextafter(problem.particleHandler.getLargestParticle()->getInteractionRadius(), std::numeric_limits<double>::max()), std::numeric_limits<double>::max());

    //for all particles, add one to the cell in cellN_ which has is associated with
    //the radius of that particle.
    for(std::vector<BaseParticle*>::const_iterator it = problem.particleHandler.begin(); it != problem.particleHandler.end(); ++it)
    {
        cellN_[radius2Cell((*it)->getInteractionRadius())]++;
    }
    
    //assign what the number of particles is.
    unsigned int numParticles = problem.particleHandler.getNumberOfObjects();

    intCellN.push_back(1.5 * cellN_[0] - 0.5 * cellN_[1]);
    for(std::vector<double>::iterator it = cellN_.begin(); it != cellN_.end(); ++it)
    {
        intCellN.push_back(*it);
    }
    intCellN.push_back(1.5 * cellN_[numCells_ - 1] - 0.5 * cellN_[numCells_ - 2]);

    if(verbosity > 0)
    {
        std::cout << "rMin=" << rMin_ << " rMax=" << rMax_ << " NParticles=" << numParticles << std::endl;
        for(unsigned int i = 0; i < cellN_.size(); i++)
        {
            std::cout << "From" << std::setw(8) << cell2Min(i) << " to" << std::setw(8) << cell2Max(i) << " number=" << std::setw(5) << cellN_[i] << std::endl;
        }
        std::cout << "Domain size [" << problem.getXMin() << "," << problem.getXMax() << "]x[" << problem.getYMin() << "," << problem.getYMax() << "]x[" << problem.getZMin() << "," << problem.getZMax() << "]" << std::endl;

        for(unsigned int i = 0; i < intCellN.size() - 1; i++)
        {
            std::cout << "From" << std::setw(8) << intCell2Min(i) << " to" << std::setw(8) << intCell2Max(i) << " linear from " << std::setw(8) << intCellN[i] << " to " << std::setw(8) << intCellN[i + 1] << std::endl;
        }
        /*std::cout<<"["<<cellN_[0];
         for (unsigned int i=1;i<cellN_.size();i++)
         {
         std::cout<<","<<cellN_[i];
         }       
         std::cout<<"]"<<std::endl;*/
    }

    dimension_ = problem.getSystemDimensions();
    switch(dimension_)
    {
    case 1:
        length_ = pow((problem.getXMax() - problem.getXMin()) / numParticles, 1);
        break;
    case 2:
        length_ = pow((problem.getXMax() - problem.getXMin()) * (problem.getYMax() - problem.getYMin()) / numParticles, 1.0 / 2.0);
        break;
    case 3:
        length_ = pow((problem.getXMax() - problem.getXMin()) * (problem.getYMax() - problem.getYMin()) * (problem.getZMax() - problem.getZMin()) / numParticles, 1.0 / 3.0);
        break;
    }
}

void HGridOptimiser::initialisePolyFunc	(	double	omega,
		std::vector< double > &	coeff,
		unsigned int	numberOfCells,
		int	verbosity
	)

Definition at line 102 of file HGridOptimiser.cc.

References cell2Max(), cell2Min(), cellCheckOverContactCheckRatio_, cellN_, dimension_, intCell2Max(), intCell2Min(), intCellN, length_, numCells_, rMax_, and rMin_.

{
    numCells_ = numberOfCells;
    cellCheckOverContactCheckRatio_ = 0.2;
    cellN_.resize(numCells_);
    for(std::vector<double>::iterator it = cellN_.begin(); it != cellN_.end(); ++it)
    {
        *it = 0;
    }

    rMin_ = nextafter(1, 0.0);
    rMax_ = nextafter(omega, std::numeric_limits<double>::max());
    for(unsigned int i = 0; i < cellN_.size(); i++)
    {
        double start = cell2Min(i);
        double end = cell2Max(i);
        for(unsigned int j = 0; j < coeff.size(); j++)
        {
            cellN_[i] += coeff[j] / (1.0 + j) * (std::pow(end, j + 1) - std::pow(start, j + 1));
        }
    }

    intCellN.push_back(1.5 * cellN_[0] - 0.5 * cellN_[1]);
    for(std::vector<double>::iterator it = cellN_.begin(); it != cellN_.end(); ++it)
    {
        intCellN.push_back(*it);
    }
    intCellN.push_back(1.5 * cellN_[numCells_ - 1] - 0.5 * cellN_[numCells_ - 2]);

    dimension_ = 1;
    length_ = 1.0;

    if(verbosity > 0)
    {
        std::cout << "rMin=" << rMin_ << " rMax=" << rMax_ << std::endl;
        for(unsigned int i = 0; i < cellN_.size(); i++)
        {
            std::cout << "From" << std::setw(8) << cell2Min(i) << " to" << std::setw(8) << cell2Max(i) << " number=" << std::setw(5) << cellN_[i] << std::endl;
        }

        for(unsigned int i = 0; i < intCellN.size() - 1; i++)
        {
            std::cout << "From" << std::setw(8) << intCell2Min(i) << " to" << std::setw(8) << intCell2Max(i) << " linear from " << std::setw(8) << intCellN[i] << " to " << std::setw(8) << intCellN[i + 1] << std::endl;
        }
    }
}

double HGridOptimiser::intCell2Max

(

unsigned int

i

)

Definition at line 180 of file HGridOptimiser.cc.

References numCells_, rMax_, and rMin_.

Referenced by diffEndExpectedCellsIntegral(), diffHExpectedCellsIntegral(), diffHExpectedCellsIntegralCellNumerator(), diffPdfInt(), diffStartExpectedCellsIntegral(), expectedCellsIntegral(), expectedCellsIntegralCellDenominator(), expectedCellsIntegralCellNumerator(), initialise(), initialisePolyFunc(), pdfInt(), and pdfIntCell().

{
    if(i == numCells_)
        return rMax_;
    else
        return rMin_ + (rMax_ - rMin_) * (i + 0.5) / numCells_;
}

double HGridOptimiser::intCell2Min

(

unsigned int

i

)

Definition at line 172 of file HGridOptimiser.cc.

References numCells_, rMax_, and rMin_.

Referenced by diffEndExpectedCellsIntegral(), diffHExpectedCellsIntegral(), diffHExpectedCellsIntegralCellNumerator(), diffPdfInt(), diffStartExpectedCellsIntegral(), expectedCellsIntegral(), expectedCellsIntegralCellDenominator(), expectedCellsIntegralCellNumerator(), initialise(), initialisePolyFunc(), pdfInt(), and pdfIntCell().

{
    if(i == 0)
        return rMin_;
    else
        return rMin_ + (rMax_ - rMin_) * (i - 0.5) / numCells_;
}

double HGridOptimiser::pdfInt	(	double	start,
		double	end,
		int	power
	)

This function calculates: int(f(r)*r^power*dr,r=start..end)/int(f(r)*dr,r=0..omega) with r=a*r+b.

Definition at line 214 of file HGridOptimiser.cc.

References intCell2Max(), intCell2Min(), numCells_, pdfIntCell(), and radius2IntCell().

Referenced by calculateDiffWork(), calculateWork(), and histNumberParticlesPerCell().

{
    //std::cout<<"pdfInit("<<start<<","<<end<<","<<p<<");"<<std::endl;
    unsigned int startCell = radius2IntCell(start);
    unsigned int endCell = radius2IntCell(end);

    double numerator = 0;
    if(startCell == endCell)
    {
        numerator = pdfIntCell(start, end, startCell, p);
    }
    else if(endCell < startCell)
    {
        numerator = 0;
    }
    else
    {
        numerator = pdfIntCell(start, intCell2Max(startCell), startCell, p);
        //std::cout<<"Teller+="<<pdfIntCell(start,intCell2Max(startCell),startCell,p)<<std::endl;
        for(unsigned int i = startCell + 1; i < endCell; i++)
        {
            numerator += pdfIntCell(intCell2Min(i), intCell2Max(i), i, p);
            //std::cout<<"Teller+="<<pdfIntCell(intCell2Min(i),intCell2Max(i),i,p)<<std::endl;
        }
        numerator += pdfIntCell(intCell2Min(endCell), end, endCell, p);
        //std::cout<<"Teller+="<<pdfIntCell(intCell2Min(endCell),end,endCell,p)<<std::endl;
    }
    double denominator = 0;
    for(unsigned int i = 0; i <= numCells_; i++)
    {
        denominator += pdfIntCell(intCell2Min(i), intCell2Max(i), i, 0);
        //std::cout<<"Noemer+="<<pdfIntCell(intCell2Min(i),intCell2Max(i),i,0)<<std::endl;
    }
    return numerator / denominator;
}

double HGridOptimiser::pdfIntCell	(	double	start,
		double	end,
		unsigned int	i,
		int	p
	)

Definition at line 204 of file HGridOptimiser.cc.

References intCell2Max(), intCell2Min(), and intCellN.

Referenced by diffPdfInt(), and pdfInt().

{
    double a = (intCellN[i + 1] - intCellN[i]) / (intCell2Max(i) - intCell2Min(i));
    double b = (intCellN[i] * intCell2Max(i) - intCellN[i + 1] * intCell2Min(i)) / (intCell2Max(i) - intCell2Min(i));
    return(a / (p + 2) * (std::pow(end, p + 2) - std::pow(start, p + 2)) + b / (p + 1) * (std::pow(end, p + 1) - std::pow(start, p + 1)));
}

unsigned int HGridOptimiser::radius2Cell

(

double

r

)

Assigns a BaseParticle of given radius to a certain cell.

Assigns a cell to a BaseParticle with the given radius. Note that the index of the cells are linear in the radius of a particle. For example, if numCells_ = 10 and we are looking in the radius range [0,10], than the cell number of the particle with radius r is the number r rounded down to an integer.

Definition at line 156 of file HGridOptimiser.cc.

References ERROR, logger, numCells_, rMax_, and rMin_.

Referenced by initialise().

{
    if(r < rMin_ || r >= rMax_)
        logger(ERROR, "The radius % is not in the range [%,%]", r, rMin_, rMax_);
    unsigned int y = static_cast<unsigned int> (floor(numCells_ * (r - rMin_) / (rMax_ - rMin_)));
    return y;
}

unsigned int HGridOptimiser::radius2IntCell

(

double

r

)

Definition at line 164 of file HGridOptimiser.cc.

References numCells_, rMax_, and rMin_.

Referenced by diffEndExpectedCellsIntegral(), diffHExpectedCellsIntegral(), diffPdfInt(), diffStartExpectedCellsIntegral(), expectedCellsIntegral(), and pdfInt().

{
    if(r < rMin_ || r >= rMax_)
        throw;
    unsigned int y = static_cast<unsigned int> (floor(numCells_ * (r - rMin_) / (rMax_ - rMin_) + 0.5));
    return y;
}

Member Data Documentation

double HGridOptimiser::cellCheckOverContactCheckRatio_

private

The ratio of the time required for a single geometric contact detection over the time required to retrieve information from a cell. This seems to be only used for checking the effort required by the HGrid, not to compute the cell sizes.

Definition at line 198 of file HGridOptimiser.h.

Referenced by calculateDiffWork(), calculateWork(), initialise(), and initialisePolyFunc().

std::vector<double> HGridOptimiser::cellN_

private

Definition at line 207 of file HGridOptimiser.h.

Referenced by initialise(), and initialisePolyFunc().

unsigned int HGridOptimiser::dimension_

private

The dimension of the system, usually 3, sometimes 2 or 1.

Definition at line 202 of file HGridOptimiser.h.

Referenced by calculateDiffWork(), calculateWork(), histNumberParticlesPerCell(), initialise(), and initialisePolyFunc().

std::vector<double> HGridOptimiser::intCellN

private

Definition at line 208 of file HGridOptimiser.h.

Referenced by diffEndExpectedCellsIntegral(), diffHExpectedCellsIntegralCellNumerator(), diffPdfInt(), diffStartExpectedCellsIntegral(), expectedCellsIntegralCellDenominator(), expectedCellsIntegralCellNumerator(), initialise(), initialisePolyFunc(), and pdfIntCell().

double HGridOptimiser::length_

private

The weighted length of the domain.

The weighted length is computed by multiplying the lengths of all directions with each other, and then taking the appropriate type of root so that the unit is the same as that of a length (square root for 2D, cube root for 3D).

Definition at line 191 of file HGridOptimiser.h.

Referenced by calculateDiffWork(), calculateWork(), histNumberParticlesPerCell(), initialise(), and initialisePolyFunc().

unsigned int HGridOptimiser::numCells_

private

Number of cells, usually called levels in the HGrid.

Definition at line 173 of file HGridOptimiser.h.

Referenced by cell2Max(), cell2Min(), diffPdfInt(), initialise(), initialisePolyFunc(), intCell2Max(), intCell2Min(), pdfInt(), radius2Cell(), and radius2IntCell().

double HGridOptimiser::rMax_

private

Radius of the largest particle, "rounded" to the smallest double that is larger than the radius of each particle.

Definition at line 183 of file HGridOptimiser.h.

Referenced by cell2Max(), cell2Min(), getOptimalDistribution(), initialise(), initialisePolyFunc(), intCell2Max(), intCell2Min(), radius2Cell(), and radius2IntCell().

double HGridOptimiser::rMin_

private

Radius of the smallest particle, "rounded" to the largest double that is smaller than the radius of each particle.

Definition at line 178 of file HGridOptimiser.h.

Referenced by cell2Max(), cell2Min(), getOptimalDistribution(), initialise(), initialisePolyFunc(), intCell2Max(), intCell2Min(), radius2Cell(), and radius2IntCell().

---

The documentation for this class was generated from the following files:

• HGridOptimiser.h
• HGridOptimiser.cc