mathsFunc Namespace Reference

Namespace for some extra maths function that are often needed
More...

Functions
Mdouble	gamma (Mdouble gamma_in)
	This is the gamma function returns the true value for the half integer value. More...
Mdouble	beta (Mdouble z, Mdouble w)
	This is the beta function, returns the approximation based on cmath's implementation of ln(gamma) More...
Mdouble	chi_squared (Mdouble x, unsigned int k)
	This is a chi_squared function return the value x and degrees of freedom k. More...
Mdouble	chi_squared_prob (Mdouble x, unsigned int k)
	This is the function which actually gives the probability back using a chi squared test. More...
Mdouble	goldenSectionSearch (Mdouble(*function)(const Mdouble), Mdouble min, Mdouble cur, Mdouble max, Mdouble endCondition, Mdouble curVal=std::numeric_limits< Mdouble >::quiet_NaN())
	This function performs a golden section search to find the location of the minimum of a function. More...
template<typename T >
int	sign (T val)
	This is a sign function, it returns -1 for negative numbers, 1 for positive numbers and 0 for 0. More...
template<typename T >
T	square (const T val)
	squares a number More...
template<typename T >
T	cubic (const T val)
	calculates the cube of a number More...
bool	isEqual (Mdouble v1, Mdouble v2, Mdouble absError)
	Compares the difference of two Mdouble with an absolute error, useful in UnitTests. More...
bool	isEqual (Vec3D v1, Vec3D v2, Mdouble absError)
	Compares the difference of two Vec3D with an absolute error, useful in UnitTests. More...
bool	isEqual (Matrix3D m1, Matrix3D m2, Mdouble absError)
	Compares the difference of two Vec3D with an absolute error, useful in UnitTests. More...
bool	isEqual (MatrixSymmetric3D m1, MatrixSymmetric3D m2, Mdouble absError)
bool	isEqual (Quaternion v1, Quaternion v2, double absError)
template<typename T >
constexpr T	factorial (const T t)
	factorial function More...
Mdouble	sin (Mdouble x)
Mdouble	cos (Mdouble x)
Mdouble	exp (Mdouble Exponent)
Mdouble	log (Mdouble Power)
template<typename T >
T	tan (T x)
Mdouble	chebyshev (Mdouble x, const Mdouble coef[], int N)
	Namespace for evaluating the zeroth modified Bessel function of the first kind, I0(x), required in StatisticsPoint.hcc. More...
Mdouble	I0_exp (Mdouble x)
Mdouble	I0 (Mdouble x)

Detailed Description

Namespace for some extra maths function that are often needed

Function Documentation

beta()

Mdouble mathsFunc::beta	(	Mdouble	z,
		Mdouble	w
	)

This is the beta function, returns the approximation based on cmath's implementation of ln(gamma)

Computes the beta function for real numbers, based on cmath's approximation of the ln of the gamma function. See https://en.wikipedia.org/wiki/Beta_function for more details

Parameters

z	first Mdouble argument to compute the beta function of
w	second Mdouble argument to compute the beta function of

Returns

the value of beta(z,w) as an Mdouble

{
    return std::exp(std::lgamma(z) + std::lgamma(w) - std::lgamma(z + w));
}

References exp().

Referenced by helpers::computeDisptFromCollisionTimeAndRestitutionCoefficientAndTangentialRestitutionCoefficientAndEffectiveMass(), SuperQuadricParticle::computeMass(), Multipole::convertMultipoleToLocal(), InfiniteWall::getFurthestPointSuperQuadric(), SuperQuadricParticle::getVolume(), LawinenBox::LawinenBox(), main(), HertzianBSHPInteractionTwoParticleElasticCollision::padeCoefficientOfRestitution(), ParhamiMcMeekingSinterSpecies::set(), SilbertPeriodic::set_study(), SuperQuadricParticle::setBoundingRadius(), SlidingFrictionSpecies::setCollisionTimeAndNormalAndTangentialRestitutionCoefficient(), SlidingFrictionSpecies::setCollisionTimeAndNormalAndTangentialRestitutionCoefficientNoDispt(), SuperQuadricParticle::setInertia(), SinterNormalSpecies::setParhamiMcKeeping(), SilbertPeriodic::SilbertPeriodic(), LocalExpansion::translateLocalExpansion(), and Multipole::TranslateMultipoleExpansionTo().

chebyshev()

Mdouble mathsFunc::chebyshev	(	Mdouble	x,
		const Mdouble	coef[],
		int	N
	)

Namespace for evaluating the zeroth modified Bessel function of the first kind, I0(x), required in StatisticsPoint.hcc.

{
    const Mdouble* p = coef;
    Mdouble b0 = *p++;
    Mdouble b1 = 0, b2;
    int i = N - 1;
    
    logger.assert_debug(i > 0, "i is greater than 0");
    do
    {
        b2 = b1;
        b1 = b0;
        b0 = x * b1 - b2 + *p++;
    } while (--i);
    
    return (0.5 * (b0 - b2));
}

References constants::i, and logger.

Referenced by I0_exp().

chi_squared()

Mdouble mathsFunc::chi_squared	(	Mdouble	x,
		unsigned int	k
	)

This is a chi_squared function return the value x and degrees of freedom k.

This is a chi_squared function return the value x and degrees of freedom k

{
    
    Mdouble prefactor = pow(2, k / 2.0) * gamma(k / 2.0);
    Mdouble mainfactor = pow(x, k / 2.0 - 1) * exp(x / -2.0);
    
    return mainfactor / prefactor;
    
}

References exp(), and gamma().

Referenced by chi_squared_prob().

chi_squared_prob()

Mdouble mathsFunc::chi_squared_prob	(	Mdouble	x,
		unsigned int	k
	)

This is the function which actually gives the probability back using a chi squared test.

This calculates the probability based on a chi squared test First we calculated the cumulative chi_squared function. This is the function which actually gives the probability back It is calculated by calling the normal chi_squared function and using the trapezoidal rule. The final results is 1-the cumulative chi_squared function

{
 
//The current value was picked by tried were it stopped effect the 4 d.p.
    const int num_steps_per_unit = 100;
    Mdouble sum = 0;
    Mdouble x = 0;
    long int num_steps = static_cast<int>(num_steps_per_unit * x_max);
//Use trapezional rule, but ignoring the ends
    for (int i = 0; i < num_steps; i++)
    {
        x = x_max / num_steps * (i + 0.5);
        sum = sum + chi_squared(x, k);
    }
    return 1.0 - sum * x_max / num_steps;
    
}

References chi_squared(), and constants::i.

Referenced by main(), and RNG::test().

cos()

Mdouble mathsFunc::cos

(

Mdouble

x

)

{
    Mdouble N = floor(x / (2.0 * constants::pi) + 0.5);
    x -= N * 2.0 * constants::pi;
    Mdouble Sum = 1.0;
    Mdouble Power = 1;
    Mdouble Sign = 1;
    const Mdouble x2 = x * x;
    Mdouble Fact = 1.0;
    for (unsigned int i = 2; i < 25; i += 2)
    {
        Power *= x2;
        Fact *= i * (i - 1);
        Sign *= -1.0;
        Sum += Sign * Power / Fact;
    }
    return Sum;
}

References constants::i, and constants::pi.

Referenced by VerticalMixerAngledBlades::addBlades(), ChuteWithHopper::addHopper(), VerticalMixerAngledBlades::addPrettyBlades(), ChuteWithWedge::addWedge(), AngledPeriodicBoundarySecondUnitTest::AngledPeriodicBoundarySecondUnitTest(), AngledPeriodicBoundaryUnitTest::AngledPeriodicBoundaryUnitTest(), LiquidBridgeWilletInteraction::computeAdhesionForce(), LiquidMigrationLSInteraction::computeAdhesionForce(), LiquidMigrationWilletInteraction::computeAdhesionForce(), BaseCluster::computeInternalStructure(), SinterLinInteraction::computeNormalForce(), Funnel::create_funnel(), SegregationWithHopper::create_inflow_particle(), AxisymmetricHopper::create_inflow_particle(), Funnel::create_inflow_particle(), CircularPeriodicBoundary::createPeriodicParticle(), LevelSetWall::createVTKSphere(), HorizontalBaseScrew::getDistanceAndNormal(), HorizontalScrew::getDistanceAndNormal(), Coil::getDistanceAndNormal(), InfiniteWallWithHole::getDistanceAndNormal(), SineWall::getDistanceAndNormal(), VChute::getDistanceAndNormal(), Screw::getDistanceAndNormalLabCoordinates(), Vec3D::getFromCylindricalCoordinates(), RNG::getNormalVariate(), vibratedBed::getPrescribedVelocity(), main(), objectivenessTest(), BaseCluster::particleInsertionSuccessful(), particleParticleTest(), HopperInsertionBoundary::placeParticle(), protectiveWall::protectiveWall(), CircularPeriodicBoundary::rotateParticle(), ClumpParticle::rotatePrincipalDirections(), ChutePeriodic::set_chute_parameters(), Cstatic3D::set_particle_properties(), Slide::set_Walls(), Quaternion::setAngleZ(), Chute::setChuteAngleAndMagnitudeOfGravity(), CurvyChute::setChuteAngleAndMagnitudeOfGravity(), Quaternion::setEuler(), T_protectiveWall::setupInitialConditions(), NozzleDemo::setupInitialConditions(), HourGlass::setupInitialConditions(), MinimalExampleDrum::setupInitialConditions(), Cstatic3D::setupInitialConditions(), RotatingDrum::setupInitialConditions(), ConstantMassFlowMaserBoundaryMixedSpeciesSelfTest::setupInitialConditions(), ConstantMassFlowMaserSelfTest::setupInitialConditions(), NozzleSelfTest::setupInitialConditions(), DrumRot::setupInitialConditions(), protectiveWall::setupInitialConditions(), Tutorial11::setupInitialConditions(), DrivenParticleClass::setupInitialConditions(), sphericalHarmonics::sphericalHarmonics(), T_protectiveWall::T_protectiveWall(), tan(), UniformRandomPDs(), Funnel::update_funnel(), wallParticleTest(), AxisymmetricIntersectionOfWalls::writeVTK(), HorizontalBaseScrew::writeVTK(), HorizontalScrew::writeVTK(), Screw::writeVTK(), and ScrewsymmetricIntersectionOfWalls::writeVTK().

cubic()

template<typename T >

T mathsFunc::cubic

(

const T

val

)

calculates the cube of a number

{
    return val * val * val;
}

Referenced by ChutePeriodic::add_flow_particles(), Vreman::add_particles(), Chute::addFlowParticlesCompactly(), CGCoordinates::XYZ::getGaussPrefactor(), PSDContinuous::getVolumetricMean(), PSD::getVolumetricMeanRadius(), LawinenBox::LawinenBox(), main(), CGCoordinates::XYZ::normalisePolynomialCoefficients(), ParticleBeam::ParticleBeam(), FlowRule::set_Polydispersity(), FlowRule::setPolydispersity(), StressStrainControl::setupInitialConditions(), AngleOfRepose::setupInitialConditions(), DPM::setupInitialConditions(), SinterForceUnitTest::setupInitialConditions(), ChuteBottom::setupInitialConditions(), ScaleCoupledBeam::setupMercury(), Siegen::Siegen(), and StressStrainControl::StressStrainControl().

exp()

Mdouble mathsFunc::exp

(

Mdouble

Exponent

)

{
    Mdouble X, P, Frac, I, L;
    X = Exponent;
    Frac = X;
    P = (1.0 + X);
    I = 1.0;
    
    do
    {
        I++;
        Frac *= (X / I);
        L = P;
        P += Frac;
    } while (L != P);
    
    return P;
}

References Global_Physical_Variables::P, and X.

Referenced by beta(), chi_squared(), HeatFluidCoupledSpecies< NormalForceSpecies >::f(), CGCoordinates::Base_X_Y_Z::getGaussIntegralPrefactor(), CGCoordinates::Base_XY_XZ_YZ::getGaussIntegralPrefactor(), CGCoordinates::XYZ::getGaussIntegralPrefactor(), CGCoordinates::Base_XY_XZ_YZ::getGaussPrefactor(), CGCoordinates::XYZ::getGaussPrefactor(), helpers::getMaximumVelocity(), RNG::getPoissonVariate(), LinearPlasticViscoelasticNormalSpecies::getRestitutionCoefficient(), LinearViscoelasticNormalSpecies::getRestitutionCoefficient(), SinterLinNormalSpecies::getRestitutionCoefficient(), SPHNormalSpecies::getRestitutionCoefficient(), HGridOptimiser::histNumberParticlesPerCell(), I0(), log(), main(), rapidjson::GenericReader< Encoding, Allocator >::ParseNumber(), particleParticleTest(), ParhamiMcMeekingSinterSpecies::set(), PSD::setDistributionLogNormal(), SinterNormalSpecies::setParhamiMcKeeping(), Calibration::setPSD(), Material::setPSD(), SilbertHstop::solve_analytic(), sphericalHarmonics::sphericalHarmonics(), and wallParticleTest().

factorial()

template<typename T >

constexpr T mathsFunc::factorial ( const T  t )

constexpr

factorial function

{
    return (t == 0) ? 1 : t * factorial(t - 1);
}

Referenced by sphericalHarmonics::computeSquaredFactorialValues(), HGridOptimiser::histNumberParticlesPerCell(), and sphericalHarmonics::sphericalHarmonics().

gamma()

Mdouble mathsFunc::gamma

(

Mdouble

gamma_in

)

This is the gamma function returns the true value for the half integer value.

This is the gamma function, gives 'exact' answers for the half integer values This is done using the recussion relation and the known values for 1 and 0.5 Note, return NaN for non-half integer values

{
    const Mdouble ep = 1e-5;
    
    if (gamma_in > 1.0 + ep)
    {
        return ((gamma_in - 1) * gamma(gamma_in - 1));
    }
    else
    {
        
        if ((gamma_in - ep < 1.0) && (gamma_in + ep > 1.0))
            return 1.0;
        else if ((gamma_in - ep < 0.5) && (gamma_in + ep > 0.5))
            return constants::sqrt_pi;
        else
            return std::numeric_limits<Mdouble>::quiet_NaN();
    }
} //end func gamma

References constants::sqrt_pi.

Referenced by chi_squared(), SegregationWithHopper::create_inflow_particle(), AngleOfRepose::create_inflow_particle(), oomph::RefineableQDPVDElement< DIM, NNODE_1D >::fill_in_generic_contribution_to_residuals_pvd(), oomph::ScaleCoupledElement< ELEMENT >::fill_in_generic_contribution_to_residuals_pvd(), oomph::VolumeCoupledElement< ELEMENT >::fill_in_generic_contribution_to_residuals_pvd(), oomph::SCoupledElement< ELEMENT >::get_momentum_and_energy(), oomph::VolumeCoupledElement< ELEMENT >::get_momentum_and_energy(), InfiniteWall::getFurthestPointSuperQuadric(), SolidProblem< ELEMENT_TYPE >::getMassMomentumEnergy(), main(), HopperInsertionBoundary::placeParticle(), and SuperQuadricParticle::setBoundingRadius().

goldenSectionSearch()

Mdouble mathsFunc::goldenSectionSearch	(	Mdouble(*)(const Mdouble)	function,
		Mdouble	min,
		Mdouble	cur,
		Mdouble	max,
		Mdouble	endCondition,
		Mdouble	curVal = std::numeric_limits<Mdouble>::quiet_NaN()
	)

This function performs a golden section search to find the location of the minimum of a function.

Parameters

[in]	function	A function pointer to the function of which you want to calculate the location of its minimum.
[in]	min	The minimum location
[in]	cur	The current location
[in]	max	The maximum location
[in]	endCondition	The algorithm terminates when abs(max - min) < endCondition
[in]	curVal	The value of the function at the current location (on default this value is calculated internally)

{
    if (std::abs(max - min) < endCondition)
    {
        return 0.5 * (min + max);
    }
    logger(INFO, "Min=% Max=% diff= %", min, max, max - min);
    Mdouble resphi = 2 - 0.5 * (1 + std::sqrt(5));
    Mdouble x;
    if (max - cur > cur - min)
    {
        x = cur + resphi * (max - cur);
    }
    else
    {
        x = cur - resphi * (cur - min);
    }
    if (std::isnan(curVal))
        curVal = function(cur);
    Mdouble xVal = function(x);
    if (xVal < curVal)
    {
        if (max - cur > cur - min)
        {
            return goldenSectionSearch(function, cur, x, max, endCondition, xVal);
        }
        else
        {
            return goldenSectionSearch(function, min, x, cur, endCondition, xVal);
        }
    }
    else
    {
        if (max - cur > cur - min)
        {
            return goldenSectionSearch(function, min, cur, x, endCondition, curVal);
        }
        else
        {
            return goldenSectionSearch(function, x, cur, max, endCondition, curVal);
        }
    }
}

References INFO, and logger.

I0()

Mdouble mathsFunc::I0

(

Mdouble

x

)

{
    if (x < 0)
        x = -x;
    return exp(x) * I0_exp(x);
}

References exp(), and I0_exp().

I0_exp()

Mdouble mathsFunc::I0_exp

(

Mdouble

x

)

{
    // Coefficients for [0..8]
    const Mdouble A[] =
            {
                    -4.415341646479339379501E-18,
                    3.330794518822238097831E-17,
                    -2.431279846547954693591E-16,
                    1.715391285555133030611E-15,
                    -1.168533287799345168081E-14,
                    7.676185498604935616881E-14,
                    -4.856446783111929460901E-13,
                    2.955052663129639834611E-12,
                    -1.726826291441555707231E-11,
                    9.675809035373236912241E-11,
                    -5.189795601635262906661E-10,
                    2.659823724682386650351E-9,
                    -1.300025009986248042121E-8,
                    6.046995022541918949321E-8,
                    -2.670793853940611733911E-7,
                    1.117387539120103718151E-6,
                    -4.416738358458750563591E-6,
                    1.644844807072889708931E-5,
                    -5.754195010082103703981E-5,
                    1.885028850958416557291E-4,
                    -5.763755745385823658851E-4,
                    1.639475616941335798421E-3,
                    -4.324309995050575944301E-3,
                    1.054646039459499831831E-2,
                    -2.373741480589946881561E-2,
                    4.930528423967070848781E-2,
                    -9.490109704804764442101E-2,
                    1.716209015222087753491E-1,
                    -3.046826723431983986831E-1,
                    6.767952744094760849951E-1
            };
    
    // Coefficients for [8..infinity]
    const Mdouble B[] =
            {
                    -7.233180487874753954561E-18,
                    -4.830504485944182071261E-18,
                    4.465621420296759999011E-17,
                    3.461222867697461093101E-17,
                    -2.827623980516583484941E-16,
                    -3.425485619677219134621E-16,
                    1.772560133056526383601E-15,
                    3.811680669352622420751E-15,
                    -9.554846698828307648701E-15,
                    -4.150569347287222086631E-14,
                    1.540086217521409826911E-14,
                    3.852778382742142701141E-13,
                    7.180124451383666233671E-13,
                    -1.794178531506806117781E-12,
                    -1.321581184044771311881E-11,
                    -3.149916527963241364541E-11,
                    1.188914710784643834241E-11,
                    4.940602388224969589101E-10,
                    3.396232025708386345151E-9,
                    2.266668990498178064591E-8,
                    2.048918589469063741831E-7,
                    2.891370520834756482971E-6,
                    6.889758346916823984261E-5,
                    3.369116478255694089901E-3,
                    8.044904110141088316081E-1
            };
    
    if (x < 0)
        x = -x;
    
    if (x <= 8.0)
    {
        Mdouble y = (x / 2.0) - 2.0;
        return (chebyshev(y, A, 30));
    }
    
    return (chebyshev(32.0 / x - 2.0, B, 25) / sqrt(x));
    
}

References A, and chebyshev().

Referenced by I0().

isEqual() [1/5]

bool mathsFunc::isEqual	(	Matrix3D	m1,
		Matrix3D	m2,
		Mdouble	absError
	)

Compares the difference of two Vec3D with an absolute error, useful in UnitTests.

Parameters

[in]	m1	The first Vec3D
[in]	m2	The second Vec3D
[in]	absError	The allowed maximum absolute error

Returns

true if the two Vec3D are equal

{
    return (isEqual(m1.XX, m2.XX, absError)
            && isEqual(m1.XY, m2.XY, absError)
            && isEqual(m1.XZ, m2.XZ, absError)
            && isEqual(m1.YX, m2.YX, absError)
            && isEqual(m1.YY, m2.YY, absError)
            && isEqual(m1.YZ, m2.YZ, absError)
            && isEqual(m1.ZX, m2.ZX, absError)
            && isEqual(m1.ZY, m2.ZY, absError)
            && isEqual(m1.ZZ, m2.ZZ, absError));
}

References isEqual(), Matrix3D::XX, Matrix3D::XY, Matrix3D::XZ, Matrix3D::YX, Matrix3D::YY, Matrix3D::YZ, Matrix3D::ZX, Matrix3D::ZY, and Matrix3D::ZZ.

isEqual() [2/5]

bool mathsFunc::isEqual	(	MatrixSymmetric3D	m1,
		MatrixSymmetric3D	m2,
		Mdouble	absError
	)

{
    return (isEqual(m1.XX, m2.XX, absError)
            && isEqual(m1.XY, m2.XY, absError)
            && isEqual(m1.XZ, m2.XZ, absError)
            && isEqual(m1.YY, m2.YY, absError)
            && isEqual(m1.YZ, m2.YZ, absError)
            && isEqual(m1.ZZ, m2.ZZ, absError));
}

References isEqual(), MatrixSymmetric3D::XX, MatrixSymmetric3D::XY, MatrixSymmetric3D::XZ, MatrixSymmetric3D::YY, MatrixSymmetric3D::YZ, and MatrixSymmetric3D::ZZ.

isEqual() [3/5]

bool mathsFunc::isEqual	(	Mdouble	v1,
		Mdouble	v2,
		Mdouble	absError
	)

Compares the difference of two Mdouble with an absolute error, useful in UnitTests.

Parameters

[in]	v1	The first Mdouble
[in]	v2	The second Mdouble
[in]	absError	The allowed maximum absolute error

Returns

True if the two Mdouble are equal

{
    return std::abs(v1 - v2) <= absError;
}

Referenced by SphericalSuperQuadricCollision::actionsAfterSolve(), ContactDetectionNormalSpheresTest::actionsAfterTimeStep(), ContactDetectionRotatedSpheresTest::actionsAfterTimeStep(), checkTemplate(), compareParticles(), isEqual(), main(), TriangleMeshWall::moveVerticesToMatchVolume(), SuperQuadricParticle::setBoundingRadius(), MercuryBase::setHGridCellOverSizeRatio(), BouncingSuperQuadric::setupInitialConditions(), SlidingSpheresUnitTest::setupInitialConditions(), ContactDetectionNormalSpheresTest::setupInitialConditions(), ContactDetectionRotatedSpheresTest::setupInitialConditions(), BoundingRadiusTester::test(), ContactDetectionNormalSpheresTest::test(), VolumeTest::test(), Packing::test(), HertzContactRestitutionUnitTest::test(), test1(), test2(), ShapeGradientHessianTester::testCushion(), ShapeGradientHessianTester::testEllipsoid(), ContactDetectionWithWallTester::testEllipsoidsContact(), ShapeGradientHessianTester::testRoundedBeam(), ShapeGradientHessianTester::testSphere(), and ContactDetectionWithWallTester::testSpheresContact().

isEqual() [4/5]

bool mathsFunc::isEqual	(	Quaternion	v1,
		Quaternion	v2,
		double	absError
	)

{
    return isEqual(v1.getComponent(0), v2.getComponent(0), absError) &&
           isEqual(v1.getComponent(1), v2.getComponent(1), absError) &&
           isEqual(v1.getComponent(2), v2.getComponent(2), absError) &&
           isEqual(v1.getComponent(3), v2.getComponent(3), absError);
}

References Quaternion::getComponent(), and isEqual().

isEqual() [5/5]

bool mathsFunc::isEqual	(	Vec3D	v1,
		Vec3D	v2,
		Mdouble	absError
	)

Compares the difference of two Vec3D with an absolute error, useful in UnitTests.

Parameters

[in]	v1	The first Vec3D
[in]	v2	The second Vec3D
[in]	absError	The allowed maximum absolute error

Returns

true if the two Vec3D are equal

{
    return isEqual(v1.X, v2.X, absError) && isEqual(v1.Y, v2.Y, absError) && isEqual(v1.Z, v2.Z, absError);
}

References isEqual(), Vec3D::X, Vec3D::Y, and Vec3D::Z.

log()

Mdouble mathsFunc::log

(

Mdouble

Power

)

Todo:

check if this function works

{
    Mdouble N, P, L, R, A, E;
    E = 2.71828182845905;
    P = Power;
    N = 0.0;
    
    // This speeds up the convergence by calculating the integral
    while (P >= E)
    {
        P /= E;
        N++;
    }
    N += (P / E);
    P = Power;
    do
    {
        A = N;
        L = (P / (exp(N - 1.0)));
        R = ((N - 1.0) * E);
        N = ((L + R) / E);
    } while (N < A);
    //} while (N != A);
    
    return N;
}

References A, Global_Physical_Variables::E, exp(), Global_Physical_Variables::P, and R.

Referenced by CircularPeriodicBoundary::checkBoundaryAfterParticleMoved(), LiquidMigrationLSInteraction::computeAdhesionForce(), helpers::computeDisptFromCollisionTimeAndRestitutionCoefficientAndTangentialRestitutionCoefficientAndEffectiveMass(), CircularPeriodicBoundary::createPeriodicParticle(), MD_demo::DistInt(), HeatFluidCoupledSpecies< NormalForceSpecies >::f(), RNG::getNormalVariate(), main(), particleParticleTest(), SlidingFrictionSpecies::setCollisionTimeAndNormalAndTangentialRestitutionCoefficient(), LinearPlasticViscoelasticNormalSpecies::setCollisionTimeAndRestitutionCoefficient(), LinearViscoelasticNormalSpecies::setCollisionTimeAndRestitutionCoefficient(), SinterLinNormalSpecies::setCollisionTimeAndRestitutionCoefficient(), SinterNormalSpecies::setCollisionTimeAndRestitutionCoefficient(), SPHNormalSpecies::setCollisionTimeAndRestitutionCoefficient(), Calibration::setPSD(), Material::setPSD(), LinearPlasticViscoelasticNormalSpecies::setRestitutionCoefficient(), LinearViscoelasticNormalSpecies::setRestitutionCoefficient(), SinterLinNormalSpecies::setRestitutionCoefficient(), SPHNormalSpecies::setRestitutionCoefficient(), SinterLinNormalSpecies::setStiffnessAndRestitutionCoefficient(), SinterNormalSpecies::setStiffnessAndRestitutionCoefficient(), StressStrainControl::setupInitialConditions(), SegregationPeriodic::setupInitialConditions(), RandomClusterInsertionBoundarySelfTest::setupInitialConditions(), and wallParticleTest().

sign()

template<typename T >

int mathsFunc::sign

(

T

val

)

This is a sign function, it returns -1 for negative numbers, 1 for positive numbers and 0 for 0.

{
    return (T(0) < val) - (val < T(0));
}

Referenced by PeriodicBoundaryHandler::checkIsReal(), SuperQuadricParticle::computeHessianLabFixed(), SuperQuadricParticle::computeShapeGradientLabFixed(), HorizontalScrew::getDistanceAndNormal(), Screw::getDistanceAndNormalLabCoordinates(), InfiniteWall::getFurthestPointSuperQuadric(), and PeriodicBoundaryHandler::updateParticles().

sin()

Mdouble mathsFunc::sin

(

Mdouble

x

)

{
    Mdouble N = floor(x / (2.0 * constants::pi) + 0.5);
    x -= N * 2.0 * constants::pi;
    Mdouble Sum = 0;
    Mdouble Power = x;
    Mdouble Sign = 1;
    const Mdouble x2 = x * x;
    Mdouble Fact = 1.0;
    for (unsigned int i = 1; i < 25; i += 2)
    {
        Sum += Sign * Power / Fact;
        Power *= x2;
        Fact *= (i + 1) * (i + 2);
        Sign *= -1.0;
    }
    return Sum;
}

References constants::i, and constants::pi.

Referenced by VerticalMixerAngledBlades::addBlades(), ChuteWithHopper::addHopper(), VerticalMixerAngledBlades::addPrettyBlades(), ChuteWithWedge::addWedge(), BaseCluster::computeInternalStructure(), SinterLinInteraction::computeNormalForce(), Funnel::create_funnel(), SegregationWithHopper::create_inflow_particle(), AxisymmetricHopper::create_inflow_particle(), Funnel::create_inflow_particle(), LevelSetWall::createVTKSphere(), HorizontalBaseScrew::getDistanceAndNormal(), HorizontalScrew::getDistanceAndNormal(), Coil::getDistanceAndNormal(), InfiniteWallWithHole::getDistanceAndNormal(), SineWall::getDistanceAndNormal(), VChute::getDistanceAndNormal(), Screw::getDistanceAndNormalLabCoordinates(), Vec3D::getFromCylindricalCoordinates(), ChuteWithHopper::getMaximumVelocityInducedByGravity(), RNG::getNormalVariate(), main(), objectivenessTest(), BaseCluster::particleInsertionSuccessful(), HopperInsertionBoundary::placeParticle(), protectiveWall::protectiveWall(), CircularPeriodicBoundary::rotateParticle(), ClumpParticle::rotatePrincipalDirections(), ScaleCoupledBeam::ScaleCoupledBeam(), ChutePeriodic::set_chute_parameters(), Cstatic3D::set_particle_properties(), Slide::set_Walls(), Quaternion::setAngleZ(), Chute::setChuteAngleAndMagnitudeOfGravity(), CurvyChute::setChuteAngleAndMagnitudeOfGravity(), Quaternion::setEuler(), T_protectiveWall::setupInitialConditions(), AxisymmetricHopper::setupInitialConditions(), NozzleDemo::setupInitialConditions(), Binary::setupInitialConditions(), my_problem::setupInitialConditions(), HourGlass::setupInitialConditions(), MinimalExampleDrum::setupInitialConditions(), Cstatic3D::setupInitialConditions(), RotatingDrum::setupInitialConditions(), ConstantMassFlowMaserBoundaryMixedSpeciesSelfTest::setupInitialConditions(), ConstantMassFlowMaserSelfTest::setupInitialConditions(), NozzleSelfTest::setupInitialConditions(), DrumRot::setupInitialConditions(), protectiveWall::setupInitialConditions(), Tutorial11::setupInitialConditions(), DrivenParticleClass::setupInitialConditions(), ExtremeOverlapWithWallsUnitTest::setupInitialConditions(), T_protectiveWall::T_protectiveWall(), tan(), UniformRandomPDs(), Funnel::update_funnel(), AxisymmetricIntersectionOfWalls::writeVTK(), HorizontalBaseScrew::writeVTK(), HorizontalScrew::writeVTK(), Screw::writeVTK(), and ScrewsymmetricIntersectionOfWalls::writeVTK().

square()

template<typename T >

T mathsFunc::square

(

const T

val

)

squares a number

{
    return val * val;
}

Referenced by AngleOfRepose::actionsBeforeTimeStep(), DPMBase::checkParticleForInteractionLocal(), helpers::computeDisptFromCollisionTimeAndRestitutionCoefficientAndTangentialRestitutionCoefficientAndEffectiveMass(), FrictionInteraction::computeFrictionForce(), MindlinRollingTorsionInteraction::computeFrictionForce(), SlidingFrictionInteraction::computeFrictionForce(), ChuteWithPeriodicInflow::computeInternalForces(), BaseParticle::computeMass(), SinterInteraction::computeNormalForce(), HertzianSinterInteraction::computeSinterForce(), LinearPlasticViscoelasticNormalSpecies::computeTimeStep(), LinearViscoelasticNormalSpecies::computeTimeStep(), SinterLinNormalSpecies::computeTimeStep(), SinterNormalSpecies::computeTimeStep(), PSD::convertProbabilityDensityNumberDistributionToProbabilityDensityVolumeDistribution(), PSD::convertProbabilityDensityToProbabilityDensityNumberDistribution(), PSDContinuous::convertSubtractiveNumberToVolume(), PSDContinuous::convertSubtractiveVolumeToNumber(), AngleOfRepose::create_inflow_particle(), Slide::create_rough_wall(), Matrix3D::deviator(), statistics_while_running< T >::finishStatistics(), SinterNormalSpecies::getCollisionTime(), LinearPlasticViscoelasticNormalSpecies::getCollisionTime(), LinearViscoelasticNormalSpecies::getCollisionTime(), SinterLinNormalSpecies::getCollisionTime(), SPHNormalSpecies::getCollisionTime(), HertzianViscoelasticNormalSpecies::getCollisionTime(), StatisticsVector< T >::getCutoff2(), HorizontalScrew::getDistanceAndNormal(), Coil::getDistanceAndNormal(), ParabolaChute::getDistanceAndNormal(), SineWall::getDistanceAndNormal(), Screw::getDistanceAndNormalLabCoordinates(), CGCoordinates::R::getDistanceSquared(), CGCoordinates::RZ::getDistanceSquared(), CGCoordinates::X::getDistanceSquared(), CGCoordinates::XY::getDistanceSquared(), CGCoordinates::XZ::getDistanceSquared(), CGCoordinates::Y::getDistanceSquared(), CGCoordinates::YZ::getDistanceSquared(), CGCoordinates::Z::getDistanceSquared(), ChargedBondedInteraction::getElasticEnergy(), IrreversibleAdhesiveInteraction::getElasticEnergy(), HertzianSinterInteraction::getElasticEnergy(), HertzianViscoelasticInteraction::getElasticEnergy(), LinearPlasticViscoelasticInteraction::getElasticEnergy(), LinearViscoelasticInteraction::getElasticEnergy(), SinterInteraction::getElasticEnergy(), SinterLinInteraction::getElasticEnergy(), SPHInteraction::getElasticEnergy(), HertzianViscoelasticInteraction::getElasticEnergyAtEquilibrium(), CGCoordinates::Base_X_Y_Z::getGaussIntegralPrefactor(), CGCoordinates::Base_XY_XZ_YZ::getGaussIntegralPrefactor(), CGCoordinates::XYZ::getGaussIntegralPrefactor(), CGCoordinates::Base_XY_XZ_YZ::getGaussPrefactor(), CGCoordinates::XYZ::getGaussPrefactor(), WallSpecies::getInfo(), helpers::getMaximumVelocity(), CGFields::GradVelocityField::getSquared(), CGFields::LiquidMigrationFields::getSquared(), CGFields::StandardFields::getSquared(), CGCoordinates::RZ::getTangentialSquared(), CGCoordinates::XY::getTangentialSquared(), CGCoordinates::XYZ::getTangentialSquared(), CGCoordinates::XZ::getTangentialSquared(), CGCoordinates::YZ::getTangentialSquared(), Mercury3D::hGridFindContactsWithTargetCell(), MercuryBase::hGridUpdateMove(), BaseParticle::isInContactWith(), main(), CGCoordinates::Base_XY_XZ_YZ::normalisePolynomialCoefficients(), particleParticleTest(), statistics_while_running< T >::processStatistics(), Vreman::set_symmetric_contraction(), ContractionWithPeriodicInflow::set_symmetric_contraction(), ChuteWithPeriodicInflowAndContraction::set_symmetric_contraction(), ChuteWithContraction::set_symmetric_contraction(), SuperQuadricParticle::setBoundingRadius(), StatisticsVector< T >::setCGWidth(), SlidingFrictionSpecies::setCollisionTimeAndNormalAndTangentialRestitutionCoefficient(), SlidingFrictionSpecies::setCollisionTimeAndNormalAndTangentialRestitutionCoefficientNoDispt(), LinearPlasticViscoelasticNormalSpecies::setCollisionTimeAndRestitutionCoefficient(), LinearViscoelasticNormalSpecies::setCollisionTimeAndRestitutionCoefficient(), SinterLinNormalSpecies::setCollisionTimeAndRestitutionCoefficient(), SinterNormalSpecies::setCollisionTimeAndRestitutionCoefficient(), SPHNormalSpecies::setCollisionTimeAndRestitutionCoefficient(), LinearPlasticViscoelasticNormalSpecies::setRestitutionCoefficient(), LinearViscoelasticNormalSpecies::setRestitutionCoefficient(), SinterLinNormalSpecies::setRestitutionCoefficient(), SPHNormalSpecies::setRestitutionCoefficient(), SinterLinNormalSpecies::setStiffnessAndRestitutionCoefficient(), SinterNormalSpecies::setStiffnessAndRestitutionCoefficient(), RotatingDrumWet::setupInitialConditions(), StressStrainControl::setupInitialConditions(), AngleOfRepose::setupInitialConditions(), HertzianSinterForceUnitTest::setupInitialConditions(), SinterForceUnitTest::setupInitialConditions(), Matrix3D::square(), MatrixSymmetric3D::square(), and wallParticleTest().

tan()

template<typename T >

T mathsFunc::tan

(

T

x

)

Todo:

should be properly computed

{
    return sin(x) / cos(x);
}

References cos(), and sin().

Referenced by ChuteWithHopper::addHopper(), BoundariesSelfTest::BoundariesSelfTest(), LiquidMigrationLSInteraction::computeAdhesionForce(), SinterLinInteraction::computeNormalForce(), SegregationWithHopper::create_inflow_particle(), FluxAndPeriodicBoundarySelfTest::FluxAndPeriodicBoundarySelfTest(), FluxBoundarySelfTest::FluxBoundarySelfTest(), VChute::getDistanceAndNormal(), SilbertHstop::IsAboveCurve(), LeesEdwardsDemo::LeesEdwardsDemo(), main(), ChuteWithHopper::setHopper(), SilbertHstop::solve_analytic(), TimeDependentPeriodicBoundary3DSelfTest::TimeDependentPeriodicBoundary3DSelfTest(), and TimeDependentPeriodicBoundaryTest::TimeDependentPeriodicBoundaryTest().