ElasticSuperQuadricUnitTest.cpp File Reference

#include <Species/HertzianViscoelasticMindlinSpecies.h>
#include <Walls/InfiniteWall.h>
#include "Mercury3D.h"
#include "Particles/SuperQuadricParticle.h"
#include "Species/LinearViscoelasticSpecies.h"
#include "Math/ExtendedMath.h"

Functions
int	main (int argc, char *argv[])

Function Documentation

main()

int main	(	int argc	,
		char *	argv[]
	)

Takes an oblique collision of elastic superellipsoids and checks whether the collision

• conserves momentum
• conserves kinetic+potential energy

{
    Mercury3D problem;
    problem.setName("ElasticSuperQuadricUnitTest");
    problem.setDomain({-1,-1,-1},{1,1,1});
    problem.setSaveCount(NEVER);
    problem.dataFile.setSaveCount(10);
    
    auto s = problem.speciesHandler.copyAndAddObject(LinearViscoelasticSpecies());
    s->setStiffness(2e5);
    s->setDensity(6.0/constants::pi);
    
    SuperQuadricParticle p;
    p.setSpecies(s);
    p.setAxes(.5,.5,.5);
    p.setExponents(.5,.9);
    p.setPosition({.1,0,.5});
    p.setVelocity({0,0,-.5});
    auto p0 = problem.particleHandler.copyAndAddObject(p);
    
    p.setPosition(-p.getPosition());
    p.setVelocity(-p.getVelocity());
    //p.rotate({0,constants::pi,0});
    auto p1 = problem.particleHandler.copyAndAddObject(p);
    
    problem.setTimeMax(1e-0);
    problem.setTimeStep(1e-4);
    
    const Vec3D momentum0 = problem.particleHandler.getMomentum();
    const Vec3D angularMomentum0 = problem.particleHandler.getAngularMomentum();
    const Mdouble kineticEnergy0 = problem.particleHandler.getKineticEnergy()
            +problem.particleHandler.getRotationalEnergy();
 
    // comment the  next to lines to turn on file output
    problem.setSuperquadricParticlesWriteVTK(true);
    problem.setFileType(FileType::NO_FILE);
    problem.solve();
 
    const Vec3D momentum1 = problem.particleHandler.getMomentum();
    const Vec3D angularMomentum1 = problem.particleHandler.getAngularMomentum();
    const Mdouble kineticEnergy1 = problem.particleHandler.getKineticEnergy()
                                   +problem.particleHandler.getRotationalEnergy();
 
    // Check conservation properties
    helpers::check(momentum0-momentum1,{0,0,0},1e-10*momentum0.getLength(),
            "Conservation of momentum");
    helpers::check(angularMomentum0-angularMomentum1,{0,0,0},1e-10*angularMomentum0.getLength(),
    "Con. of angular momentum");
    helpers::check(kineticEnergy0-kineticEnergy1,0,1e-4*kineticEnergy0,
            "Con. of kinetic energy  ");
    // Note, energy conservation is not very good; it does not decay quickly with time step
 
    return 0;
}

References helpers::check(), e(), Vec3D::getLength(), NEVER, NO_FILE, p, p0, p1, constants::pi, problem, and s.