CGExactOverlapUnitTestOverlap.cpp File Reference

#include "Mercury3D.h"
#include "Species/LinearViscoelasticSpecies.h"
#include "CG/TimeAveragedCG.h"
#include "CG/Functions/ExactOverlap.h"
#include "CG/Fields/StandardFieldsBinning.h"

Classes
class	CGExactOverlapUnitTest
	Tests if the displacement momentum density is correctly calculated. More...

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

Function Documentation

main()

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

{
    //declare the DPM problem and set the name
    CGExactOverlapUnitTest problem;
    problem.setTimeMax(0.1);
    problem.setTimeStep(0.1);
    problem.setFileType(FileType::NO_FILE);
 
    // Test a particle sitting at the corner of the grid cell
    Mdouble radius = 0.1;
    problem.setParameters(radius, radius*Vec3D(-1, -1, -1), radius*Vec3D(1, 1, 1), radius*Vec3D(1, 1, 1));
    auto cg0 = problem.cgHandler.copyAndAddObject(TimeAveragedCG<CGCoordinates::XYZ,CGFunctions::ExactOverlap,CGFields::StandardFieldsBinning>());
    cg0->selectSpecies(0);
    
    //run the simulation
    problem.solve(argc, argv);
 
    Mdouble density = cg0->evaluateAverage().getDensity();
    Mdouble expectedDensity = constants::pi * radius * radius * radius / 6 / problem.getVolume();
 
    logger.assert_always(std::abs(density-expectedDensity) < 1e-10, "Density is %, but should be %", density, expectedDensity);
 
    // Place a particle, so that a cap of height h is inside the grid cell
    Mdouble h = 0.05;
    radius = 0.1;
    problem.clear();
    problem.setParameters(radius, Vec3D(-2*radius, -2*radius, radius-h), Vec3D(2*radius, 2*radius, 2*radius), radius*Vec3D(0, 0, 0));
    auto cg1 = problem.cgHandler.copyAndAddObject(TimeAveragedCG<CGCoordinates::XYZ,CGFunctions::ExactOverlap,CGFields::StandardFieldsBinning>());
    cg1->selectSpecies(0);
 
    //run the simulation
    problem.solve(argc, argv);
 
    density = cg1->evaluateAverage().getDensity();
    expectedDensity = constants::pi * h * h / 3 * (3*radius - h) / problem.getVolume();
 
    logger.assert_always(std::abs(density-expectedDensity) < 1e-10, "Spherecap: Density is %, but should be %", density, expectedDensity);
 
    // Place a particle, so that a spherical segemtn with height h and distance d between the segments bottom and the sphere's center is inside the grid cell 
    h = 0.05;
    Mdouble d = -0.025;
    radius = 0.1;
 
    // The bottom radius
    Mdouble a1 = std::sqrt(radius*radius - d*d);
    // The top radius
    Mdouble a2 = std::sqrt(radius*radius - (d+h)*(d+h));
 
    logger(INFO, "a1 = %, a2 = %", a1, a2);
 
    problem.clear();
    problem.setParameters(radius, Vec3D(-2*radius, -2*radius, d), Vec3D(2*radius, 2*radius, d+h), Vec3D(0, 0, 0));
    auto cg2 = problem.cgHandler.copyAndAddObject(TimeAveragedCG<CGCoordinates::XYZ,CGFunctions::ExactOverlap,CGFields::StandardFieldsBinning>());
    cg2->selectSpecies(0);
 
    //run the simulation
    problem.solve(argc, argv);
 
    density = cg2->evaluateAverage().getDensity();
    expectedDensity = constants::pi * h * (3*a1*a1 + 3*a2*a2 + h*h) / 6 / problem.getVolume();
 
    logger.assert_always(std::abs(density-expectedDensity) < 1e-10, "Spherical segment: Density is %, but should be %", density, expectedDensity);
 
    // Place a particle, so that a spherical segemtn with height h and distance d between the segments bottom and the sphere's center is inside the grid cell 
    h = 0.1;
    d = -0.025;
    radius = 0.1;
 
    // The bottom radius
    a1 = std::sqrt(radius*radius - d*d);
    // The top radius
    a2 = std::sqrt(radius*radius - (d+h)*(d+h));
    logger(INFO, "a1 = %, a2 = %", a1, a2);
 
    problem.clear();
    problem.setParameters(radius, Vec3D(-2*radius, -2*radius, d), Vec3D(2*radius, 2*radius, d+h), Vec3D(0, 0, 0));
    auto cg3 = problem.cgHandler.copyAndAddObject(TimeAveragedCG<CGCoordinates::XYZ,CGFunctions::ExactOverlap,CGFields::StandardFieldsBinning>());
    cg3->selectSpecies(0);
 
    //run the simulation
    problem.solve(argc, argv);
 
    density = cg3->evaluateAverage().getDensity();
    expectedDensity = constants::pi * h * (3*a1*a1 + 3*a2*a2 + h*h) / 6 / problem.getVolume();
 
    logger.assert_always(std::abs(density-expectedDensity) < 1e-10, "Spherical segment 2: Density is %, but should be %", density, expectedDensity);
}

References abs(), UniformPSDSelfTest::density, e(), INFO, logger, NO_FILE, constants::pi, problem, UniformPSDSelfTest::radius, and sqrt().