QuaternionUnitTest.cpp File Reference

#include "Species/LinearViscoelasticSpecies.h"
#include "DPMBase.h"

Functions
void	test1 ()
void	test2 ()
int	main (int argc, char **argv)

Function Documentation

main()

int main	(	int argc	,
		char **	argv
	)

                                {
    test1();
    test2();
    return 0;
}

References test1(), and test2().

test1()

void test1

(

)

Todo:

shouldn't this be default?

             {
    logger(INFO, "Test 1:\n"
                 "We apply a torque T around the z-axis to a motionless\n"
                 "particle of inertia I for a time t using a time step dt.\n"
                 "The result is a final angular velocity omega=T/I*t\n"
                 "and rotation angle alpha=T/I*t^2/2.\n"
                 "\n"
                 "The input/output values should be:\n"
                 "   T = 0.1*pi*(0 0 1) Nm\n"
                 "   I = 0.1 kg m^2\n"
                 "   t = 1 s\n"
                 "   dt = 1e-5 s\n"
                 "   omega=pi*(0 0 1) rad/s\n"
                 "   alpha = pi/2*(0 0 1) rad\n");
    
    //compute mass (requires species and handler to be set)
    DPMBase D;
    LinearViscoelasticSpecies* S = D.speciesHandler.copyAndAddObject(LinearViscoelasticSpecies());
    BaseParticle* P = D.particleHandler.copyAndAddObject(SphericalParticle(S));
    D.setDimension(3); 
    S->setDensity(6.0 / constants::pi);
    P->setRadius(0.5);
    S->computeMass(P);
    
    //std::cout << "P " << *P << std::endl;
    P->setForce({1, 0, 0});
    P->setTorque(0.1 * constants::pi * Vec3D(0, 0, 1));
    Mdouble timeMax = 1.0;
    Mdouble timeStep = 1e-5;
    for (Mdouble time = 0; time < timeMax; time += timeStep)
    {
        P->accelerate(P->getForce() * P->getInvMass() * timeStep);
        P->move(P->getVelocity() * timeStep);
        P->angularAccelerate(
                P->getOrientation().rotateInverseInertiaTensor(P->getInvInertia()) * P->getTorque() * timeStep);
        P->rotate(P->getAngularVelocity() * timeStep);
    }
    
    logger(INFO, "Results:\n"
                 "   omega= %\n"
                 "   alpha= %", P->getAngularVelocity(), P->getOrientation().getEuler());
    
    //    std::cout << P->getAngularVelocity() - Vec3D(0,0,pi) << std::endl;
    //    std::cout << P->getOrientation().getEuler() - Vec3D(0,0,pi/2) << std::endl;
    
    //check results (with numerical error values added))
    if (!mathsFunc::isEqual(P->getAngularVelocity(), Vec3D(0, 0, pi + 3.14159e-05), 1e-10))
    {
        logger(ERROR, "angular velocity is %, but should be %", P->getAngularVelocity(), pi / 2);
    }
    if (!mathsFunc::isEqual(P->getOrientation().getEuler(), {0, 0, pi / 2 + 4.71241e-05}, 1e-10))
    {
        logger(ERROR, "orientation is %, but should be %", P->getOrientation().getEuler().Z, pi / 2);
    }
}

References ParticleSpecies::computeMass(), D, ERROR, INFO, mathsFunc::isEqual(), logger, Global_Physical_Variables::P, constants::pi, and ParticleSpecies::setDensity().

Referenced by main().

test2()

void test2

(

)

Todo:

shouldn't this be default?

{
    logger(INFO, "Test 2:\n"
                 "We apply a torque T around the z-axis to a motionless\n"
                 "particle of inertia I for a time t using a time step dt.\n"
                 "The result is a final angular velocity omega=inv(I)*T*t\n"
                 "and rotation angle alpha=inv(I)*T*t^2/2.\n"
                 "\n"
                 "The input/output values should be:\n"
                 "   T = 0.1*pi*(0 0 1) Nm\n"
                 "   I = 0.1*(2 0 0;0 2 0;0 0 1) kg m^2\n"
                 "   t = 1 s\n"
                 "   dt = 1e-5 s\n"
                 "   omega=pi*(0 0 1) rad/s\n"
                 "   alpha = pi/2*(0 0 1) rad\n");
    
    //compute mass (requires species and handler to be set)
    DPMBase D;
    LinearViscoelasticSpecies* S = D.speciesHandler.copyAndAddObject(LinearViscoelasticSpecies());
    BaseParticle* P = D.particleHandler.copyAndAddObject(SphericalParticle(S));
    D.setDimension(3); 
    S->setDensity(6.0 / constants::pi);
    P->setRadius(0.5);
    S->computeMass(P);
    P->setInertia(MatrixSymmetric3D(2, 1, 0,
                                    2, 0,
                                    1) * 0.1);
    
    //std::cout << "P " << *P << std::endl;
    P->setForce({1, 0, 0});
    P->setTorque(0.1 * constants::pi * Vec3D(0, 0, 1));
    Mdouble timeMax = 1.0;
    Mdouble timeStep = 1e-5;
    for (Mdouble time = 0; time < timeMax; time += timeStep)
    {
        P->accelerate(P->getForce() * P->getInvMass() * timeStep);
        P->move(P->getVelocity() * timeStep);
        P->angularAccelerate(
                P->getOrientation().rotateInverseInertiaTensor(P->getInvInertia()) * P->getTorque() * timeStep);
        P->rotate(P->getAngularVelocity() * timeStep);
    }
    
    logger(INFO, "Results:\n"
                 "   omega= %\n"
                 "   alpha= %", P->getAngularVelocity(), P->getOrientation().getEuler());
    
    //    std::cout << P->getAngularVelocity() - Vec3D(0,0,pi) << std::endl;
    //    std::cout << P->getOrientation().getEuler() - Vec3D(0,0,pi/2) << std::endl;
    
    //check results (with numerical error values added))
    if (!mathsFunc::isEqual(P->getAngularVelocity(), Vec3D(0, 0, pi + 3.14159e-05), 1e-10))
    {
        logger(ERROR, "angular velocity is %, but should be %", P->getAngularVelocity(), pi / 2);
    }
    if (!mathsFunc::isEqual(P->getOrientation().getEuler(), {0, 0, pi / 2 + 4.71241e-05}, 1e-10))
    {
        logger(ERROR, "orientation is %, but should be %", P->getOrientation().getEuler().Z, pi / 2);
    }
}

References ParticleSpecies::computeMass(), D, ERROR, INFO, mathsFunc::isEqual(), logger, Global_Physical_Variables::P, constants::pi, and ParticleSpecies::setDensity().

Referenced by main().