Commit 30a0eccf authored by Sebastian Eibl's avatar Sebastian Eibl

Merge branch 'remove_with_shape' into 'master'

[API] merge integrators with/without shape

See merge request !262
parents 32050e1e a283e38f
Pipeline #22950 failed with stages
in 680 minutes and 53 seconds
......@@ -26,8 +26,8 @@
#include "mesa_pd/data/ShapeStorage.h"
#include "mesa_pd/kernel/DoubleCast.h"
#include "mesa_pd/kernel/ExplicitEulerWithShape.h"
#include "mesa_pd/kernel/VelocityVerletWithShape.h"
#include "mesa_pd/kernel/ExplicitEuler.h"
#include "mesa_pd/kernel/VelocityVerlet.h"
#include "mesa_pd/kernel/LinearSpringDashpot.h"
#include "mesa_pd/mpi/ReduceContactHistory.h"
......@@ -130,11 +130,11 @@ int main( int argc, char ** argv )
data::particle_flags::set(p0.getFlagsRef(), data::particle_flags::FIXED);
// velocity verlet
kernel::VelocityVerletWithShapePreForceUpdate vvPreForce( dt );
kernel::VelocityVerletWithShapePostForceUpdate vvPostForce( dt );
kernel::VelocityVerletPreForceUpdate vvPreForce( dt );
kernel::VelocityVerletPostForceUpdate vvPostForce( dt );
// explicit euler
kernel::ExplicitEulerWithShape explEuler( dt );
kernel::ExplicitEuler explEuler( dt );
// collision response
collision_detection::AnalyticContactDetection acd;
......
......@@ -80,10 +80,10 @@
#include "mesa_pd/domain/BlockForestDomain.h"
#include "mesa_pd/domain/BlockForestDataHandling.h"
#include "mesa_pd/kernel/DoubleCast.h"
#include "mesa_pd/kernel/ExplicitEulerWithShape.h"
#include "mesa_pd/kernel/ExplicitEuler.h"
#include "mesa_pd/kernel/LinearSpringDashpot.h"
#include "mesa_pd/kernel/ParticleSelector.h"
#include "mesa_pd/kernel/VelocityVerletWithShape.h"
#include "mesa_pd/kernel/VelocityVerlet.h"
#include "mesa_pd/mpi/SyncNextNeighbors.h"
#include "mesa_pd/mpi/ReduceProperty.h"
#include "mesa_pd/mpi/ReduceContactHistory.h"
......@@ -853,9 +853,9 @@ int main( int argc, char **argv )
syncCall();
real_t timeStepSizeRPD = real_t(1)/real_t(numRPDSubCycles);
mesa_pd::kernel::ExplicitEulerWithShape explEulerIntegrator(timeStepSizeRPD);
mesa_pd::kernel::VelocityVerletWithShapePreForceUpdate vvIntegratorPreForce(timeStepSizeRPD);
mesa_pd::kernel::VelocityVerletWithShapePostForceUpdate vvIntegratorPostForce(timeStepSizeRPD);
mesa_pd::kernel::ExplicitEuler explEulerIntegrator(timeStepSizeRPD);
mesa_pd::kernel::VelocityVerletPreForceUpdate vvIntegratorPreForce(timeStepSizeRPD);
mesa_pd::kernel::VelocityVerletPostForceUpdate vvIntegratorPostForce(timeStepSizeRPD);
// linear model
mesa_pd::kernel::LinearSpringDashpot linearCollisionResponse(1);
......
......@@ -72,9 +72,9 @@
#include "mesa_pd/data/shape/Sphere.h"
#include "mesa_pd/domain/BlockForestDomain.h"
#include "mesa_pd/kernel/DoubleCast.h"
#include "mesa_pd/kernel/ExplicitEulerWithShape.h"
#include "mesa_pd/kernel/ExplicitEuler.h"
#include "mesa_pd/kernel/ParticleSelector.h"
#include "mesa_pd/kernel/VelocityVerletWithShape.h"
#include "mesa_pd/kernel/VelocityVerlet.h"
#include "mesa_pd/mpi/SyncNextNeighbors.h"
#include "mesa_pd/mpi/ReduceProperty.h"
#include "mesa_pd/mpi/ContactFilter.h"
......@@ -639,9 +639,9 @@ int main( int argc, char **argv )
syncCall();
mesa_pd::kernel::ExplicitEulerWithShape explEulerIntegrator(real_t(1)/real_t(numRPDSubCycles));
mesa_pd::kernel::VelocityVerletWithShapePreForceUpdate vvIntegratorPreForce(real_t(1)/real_t(numRPDSubCycles));
mesa_pd::kernel::VelocityVerletWithShapePostForceUpdate vvIntegratorPostForce(real_t(1)/real_t(numRPDSubCycles));
mesa_pd::kernel::ExplicitEuler explEulerIntegrator(real_t(1)/real_t(numRPDSubCycles));
mesa_pd::kernel::VelocityVerletPreForceUpdate vvIntegratorPreForce(real_t(1)/real_t(numRPDSubCycles));
mesa_pd::kernel::VelocityVerletPostForceUpdate vvIntegratorPostForce(real_t(1)/real_t(numRPDSubCycles));
mesa_pd::mpi::ReduceProperty reduceProperty;
......
......@@ -76,7 +76,7 @@
#include "mesa_pd/data/shape/Sphere.h"
#include "mesa_pd/domain/BlockForestDomain.h"
#include "mesa_pd/kernel/DoubleCast.h"
#include "mesa_pd/kernel/ExplicitEulerWithShape.h"
#include "mesa_pd/kernel/ExplicitEuler.h"
#include "mesa_pd/kernel/ParticleSelector.h"
#include "mesa_pd/mpi/SyncNextNeighbors.h"
#include "mesa_pd/mpi/ReduceProperty.h"
......@@ -597,7 +597,7 @@ int main( int argc, char **argv )
syncCall();
mesa_pd::kernel::ExplicitEulerWithShape explEulerIntegrator(real_t(1)/real_t(numRPDSubCycles));
mesa_pd::kernel::ExplicitEuler explEulerIntegrator(real_t(1)/real_t(numRPDSubCycles));
mesa_pd::mpi::ReduceProperty reduceProperty;
......
......@@ -80,11 +80,11 @@
#include "mesa_pd/domain/BlockForestDomain.h"
#include "mesa_pd/domain/BlockForestDataHandling.h"
#include "mesa_pd/kernel/DoubleCast.h"
#include "mesa_pd/kernel/ExplicitEulerWithShape.h"
#include "mesa_pd/kernel/ExplicitEuler.h"
#include "mesa_pd/kernel/LinearSpringDashpot.h"
#include "mesa_pd/kernel/NonLinearSpringDashpot.h"
#include "mesa_pd/kernel/ParticleSelector.h"
#include "mesa_pd/kernel/VelocityVerletWithShape.h"
#include "mesa_pd/kernel/VelocityVerlet.h"
#include "mesa_pd/mpi/ClearNextNeighborSync.h"
#include "mesa_pd/mpi/SyncNextNeighbors.h"
#include "mesa_pd/mpi/ReduceProperty.h"
......@@ -962,9 +962,9 @@ int main( int argc, char **argv )
syncCall();
real_t timeStepSizeRPD = real_t(1)/real_t(numRPDSubCycles);
mesa_pd::kernel::ExplicitEulerWithShape explEulerIntegrator(timeStepSizeRPD);
mesa_pd::kernel::VelocityVerletWithShapePreForceUpdate vvIntegratorPreForce(timeStepSizeRPD);
mesa_pd::kernel::VelocityVerletWithShapePostForceUpdate vvIntegratorPostForce(timeStepSizeRPD);
mesa_pd::kernel::ExplicitEuler explEulerIntegrator(timeStepSizeRPD);
mesa_pd::kernel::VelocityVerletPreForceUpdate vvIntegratorPreForce(timeStepSizeRPD);
mesa_pd::kernel::VelocityVerletPostForceUpdate vvIntegratorPostForce(timeStepSizeRPD);
// linear model
mesa_pd::kernel::LinearSpringDashpot linearCollisionResponse(1);
......
......@@ -50,7 +50,6 @@ if __name__ == '__main__':
mpd.add(kernel.DoubleCast(ps))
mpd.add(kernel.ExplicitEuler())
mpd.add(kernel.ExplicitEulerWithShape())
mpd.add(kernel.ForceLJ())
mpd.add(kernel.HeatConduction())
mpd.add(kernel.InsertParticleIntoLinkedCells())
......@@ -60,7 +59,7 @@ if __name__ == '__main__':
mpd.add(kernel.SpringDashpot())
mpd.add(kernel.TemperatureIntegration())
mpd.add(kernel.VelocityVerlet())
mpd.add(kernel.VelocityVerletWithShape())
mpd.add(kernel.VelocityVerlet())
mpd.add(mpi.BroadcastProperty())
mpd.add(mpi.ClearGhostOwnerSync())
......
......@@ -38,7 +38,7 @@
#include <mesa_pd/domain/BlockForestDomain.h>
#include <mesa_pd/kernel/AssocToBlock.h>
#include <mesa_pd/kernel/DoubleCast.h>
#include <mesa_pd/kernel/ExplicitEulerWithShape.h>
#include <mesa_pd/kernel/ExplicitEuler.h>
#include <mesa_pd/kernel/InsertParticleIntoLinkedCells.h>
#include <mesa_pd/kernel/ParticleSelector.h>
#include <mesa_pd/kernel/SpringDashpot.h>
......@@ -246,7 +246,7 @@ int main( int argc, char ** argv )
WALBERLA_LOG_INFO_ON_ROOT("*** SIMULATION - START ***");
// Init kernels
kernel::ExplicitEulerWithShape explicitEulerWithShape( params.dt );
kernel::ExplicitEuler explicitEuler( params.dt );
DEM dem(domain, params.dt, ss->shapes[smallSphere]->getMass());
kernel::InsertParticleIntoLinkedCells ipilc;
kernel::AssocToBlock assoc(forest);
......@@ -316,7 +316,7 @@ int main( int argc, char ** argv )
tp["Euler"].start();
//ps->forEachParticle(false, [&](const size_t idx){WALBERLA_CHECK_EQUAL(ps->getForce(idx), Vec3(0,0,0), *(*ps)[idx] << "\n" << idx);});
ps->forEachParticle(true, kernel::SelectLocal(), accessor, explicitEulerWithShape, accessor);
ps->forEachParticle(true, kernel::SelectLocal(), accessor, explicitEuler, accessor);
if (params.bBarrier) WALBERLA_MPI_BARRIER();
tp["Euler"].end();
......
......@@ -38,7 +38,7 @@
#include <mesa_pd/domain/BlockForestDomain.h>
#include <mesa_pd/kernel/AssocToBlock.h>
#include <mesa_pd/kernel/DoubleCast.h>
#include <mesa_pd/kernel/ExplicitEulerWithShape.h>
#include <mesa_pd/kernel/ExplicitEuler.h>
#include <mesa_pd/kernel/InsertParticleIntoLinkedCells.h>
#include <mesa_pd/kernel/ParticleSelector.h>
#include <mesa_pd/kernel/SpringDashpot.h>
......@@ -263,7 +263,7 @@ int main( int argc, char ** argv )
WALBERLA_LOG_INFO_ON_ROOT("*** SIMULATION - START ***");
// Init kernels
kernel::ExplicitEulerWithShape explicitEulerWithShape( params.dt );
kernel::ExplicitEuler explicitEuler( params.dt );
kernel::AssocToBlock assoc(forest);
kernel::InsertParticleIntoLinkedCells ipilc;
kernel::SpringDashpot sd(1);
......@@ -398,7 +398,7 @@ int main( int argc, char ** argv )
tp["Euler"].start();
for (int64_t i=0; i < params.simulationSteps; ++i)
{
ps->forEachParticle(true, kernel::SelectLocal(), accessor, explicitEulerWithShape, accessor);
ps->forEachParticle(true, kernel::SelectLocal(), accessor, explicitEuler, accessor);
}
tp["Euler"].end();
LIKWID_MARKER_STOP("Euler");
......
......@@ -41,7 +41,7 @@
#include <mesa_pd/domain/InfoCollection.h>
#include <mesa_pd/kernel/AssocToBlock.h>
#include <mesa_pd/kernel/DoubleCast.h>
#include <mesa_pd/kernel/ExplicitEulerWithShape.h>
#include <mesa_pd/kernel/ExplicitEuler.h>
#include <mesa_pd/kernel/InsertParticleIntoLinkedCells.h>
#include <mesa_pd/kernel/InsertParticleIntoSparseLinkedCells.h>
#include <mesa_pd/kernel/ParticleSelector.h>
......@@ -236,7 +236,7 @@ int main( int argc, char ** argv )
WALBERLA_LOG_INFO_ON_ROOT("*** SIMULATION - START ***");
// Init kernels
kernel::ExplicitEulerWithShape explicitEulerWithShape( params.dt );
kernel::ExplicitEuler explicitEuler( params.dt );
kernel::InsertParticleIntoLinkedCells ipilc;
kernel::SpringDashpot dem(1);
dem.setStiffness(0, 0, real_t(0));
......@@ -334,7 +334,7 @@ int main( int argc, char ** argv )
tpImbalanced["Euler"].start();
for (int64_t i=0; i < params.simulationSteps; ++i)
{
ps->forEachParticle(true, kernel::SelectLocal(), accessor, explicitEulerWithShape, accessor);
ps->forEachParticle(true, kernel::SelectLocal(), accessor, explicitEuler, accessor);
}
tpImbalanced["Euler"].end();
......@@ -479,7 +479,7 @@ int main( int argc, char ** argv )
tpBalanced["Euler"].start();
for (int64_t i=0; i < params.simulationSteps; ++i)
{
ps->forEachParticle(true, kernel::SelectLocal(), accessor, explicitEulerWithShape, accessor);
ps->forEachParticle(true, kernel::SelectLocal(), accessor, explicitEuler, accessor);
}
tpBalanced["Euler"].end();
......
......@@ -41,7 +41,7 @@
#include <mesa_pd/domain/InfoCollection.h>
#include <mesa_pd/kernel/AssocToBlock.h>
#include <mesa_pd/kernel/DoubleCast.h>
#include <mesa_pd/kernel/ExplicitEulerWithShape.h>
#include <mesa_pd/kernel/ExplicitEuler.h>
#include <mesa_pd/kernel/InsertParticleIntoLinkedCells.h>
#include <mesa_pd/kernel/ParticleSelector.h>
#include <mesa_pd/kernel/SpringDashpot.h>
......@@ -230,7 +230,7 @@ int main( int argc, char ** argv )
WALBERLA_LOG_INFO_ON_ROOT("*** SIMULATION - START ***");
// Init kernels
kernel::ExplicitEulerWithShape explicitEulerWithShape( params.dt );
kernel::ExplicitEuler explicitEuler( params.dt );
kernel::InsertParticleIntoLinkedCells ipilc;
kernel::SpringDashpot dem(1);
dem.setStiffness(0, 0, real_t(0));
......@@ -320,7 +320,7 @@ int main( int argc, char ** argv )
tpImbalanced["Euler"].start();
//ps->forEachParticle(false, [&](const size_t idx){WALBERLA_CHECK_EQUAL(ps->getForce(idx), Vec3(0,0,0), *(*ps)[idx] << "\n" << idx);});
ps->forEachParticle(true, kernel::SelectLocal(), accessor, explicitEulerWithShape, accessor);
ps->forEachParticle(true, kernel::SelectLocal(), accessor, explicitEuler, accessor);
if (params.bBarrier) WALBERLA_MPI_BARRIER();
tpImbalanced["Euler"].end();
......@@ -418,7 +418,7 @@ int main( int argc, char ** argv )
tpBalanced["Euler"].start();
//ps->forEachParticle(false, [&](const size_t idx){WALBERLA_CHECK_EQUAL(ps->getForce(idx), Vec3(0,0,0), *(*ps)[idx] << "\n" << idx);});
ps->forEachParticle(true, kernel::SelectLocal(), accessor, explicitEulerWithShape, accessor);
ps->forEachParticle(true, kernel::SelectLocal(), accessor, explicitEuler, accessor);
if (params.bBarrier) WALBERLA_MPI_BARRIER();
tpBalanced["Euler"].end();
......
......@@ -36,7 +36,7 @@
#include <mesa_pd/data/ShapeStorage.h>
#include <mesa_pd/domain/BlockForestDomain.h>
#include <mesa_pd/kernel/DoubleCast.h>
#include <mesa_pd/kernel/ExplicitEulerWithShape.h>
#include <mesa_pd/kernel/ExplicitEuler.h>
#include <mesa_pd/kernel/InsertParticleIntoLinkedCells.h>
#include <mesa_pd/kernel/ParticleSelector.h>
#include <mesa_pd/kernel/SpringDashpot.h>
......@@ -173,7 +173,7 @@ int main( int argc, char ** argv )
WALBERLA_LOG_INFO_ON_ROOT("*** KERNELS ***");
//! [Kernels]
kernel::ExplicitEulerWithShape explicitEulerWithShape( dt );
kernel::ExplicitEuler explicitEuler( dt );
kernel::InsertParticleIntoLinkedCells ipilc;
kernel::SpringDashpot dem( 2 );
auto meffSphereSphere = real_t(0.5) * (real_c(4.0)/real_c(3.0) * math::pi) * radius * radius * radius * density;
......@@ -218,7 +218,7 @@ int main( int argc, char ** argv )
RP.operator()<ForceTorqueNotification>(*ps);
ps->forEachParticle(true, kernel::SelectLocal(), accessor, explicitEulerWithShape, accessor);
ps->forEachParticle(true, kernel::SelectLocal(), accessor, explicitEuler, accessor);
SNN(*ps, domain);
}
......
......@@ -85,7 +85,6 @@ if __name__ == '__main__':
mpd.add(kernel.DetectAndStoreContacts())
mpd.add(kernel.DoubleCast(ps))
mpd.add(kernel.ExplicitEuler())
mpd.add(kernel.ExplicitEulerWithShape())
mpd.add(kernel.ForceLJ())
mpd.add(kernel.HCSITSRelaxationStep())
mpd.add(kernel.HeatConduction())
......@@ -101,7 +100,6 @@ if __name__ == '__main__':
mpd.add(kernel.SpringDashpotSpring())
mpd.add(kernel.TemperatureIntegration())
mpd.add(kernel.VelocityVerlet())
mpd.add(kernel.VelocityVerletWithShape())
mpd.add(mpi.BroadcastProperty())
mpd.add(mpi.ClearGhostOwnerSync())
......
......@@ -5,13 +5,19 @@ from mesa_pd.utility import generate_file
class ExplicitEuler:
def __init__(self):
self.context = dict()
self.context['interface'] = [create_access("position", "walberla::mesa_pd::Vec3", access="gs"),
create_access("linearVelocity", "walberla::mesa_pd::Vec3", access="gs"),
create_access("invMass", "walberla::real_t", access="g"),
create_access("force", "walberla::mesa_pd::Vec3", access="gs"),
create_access("flags", "walberla::mesa_pd::data::particle_flags::FlagT", access="g")]
def __init__(self, integrate_rotation=True):
self.context = {'bIntegrateRotation': integrate_rotation, 'interface': []}
self.context['interface'].append(create_access("position", "walberla::mesa_pd::Vec3", access="gs"))
self.context['interface'].append(create_access("linearVelocity", "walberla::mesa_pd::Vec3", access="gs"))
self.context['interface'].append(create_access("invMass", "walberla::real_t", access="g"))
self.context['interface'].append(create_access("force", "walberla::mesa_pd::Vec3", access="gs"))
self.context['interface'].append(create_access("flags", "walberla::mesa_pd::data::particle_flags::FlagT", access="g"))
if integrate_rotation:
self.context['interface'].append(create_access("rotation", "walberla::mesa_pd::Rot3", access="gs"))
self.context['interface'].append(create_access("angularVelocity", "walberla::mesa_pd::Vec3", access="gs"))
self.context['interface'].append(create_access("invInertiaBF", "walberla::mesa_pd::Mat3", access="g"))
self.context['interface'].append(create_access("torque", "walberla::mesa_pd::Vec3", access="gs"))
def generate(self, module):
ctx = {'module': module, **self.context}
......
# -*- coding: utf-8 -*-
from mesa_pd.accessor import create_access
from mesa_pd.utility import generate_file
class ExplicitEulerWithShape:
def __init__(self):
self.context = dict()
self.context['interface'] = [create_access("position", "walberla::mesa_pd::Vec3", access="gs"),
create_access("linearVelocity", "walberla::mesa_pd::Vec3", access="gs"),
create_access("invMass", "walberla::real_t", access="g"),
create_access("force", "walberla::mesa_pd::Vec3", access="gs"),
create_access("rotation", "walberla::mesa_pd::Rot3", access="gs"),
create_access("angularVelocity", "walberla::mesa_pd::Vec3", access="gs"),
create_access("invInertiaBF", "walberla::mesa_pd::Mat3", access="g"),
create_access("torque", "walberla::mesa_pd::Vec3", access="gs"),
create_access("flags", "walberla::mesa_pd::data::particle_flags::FlagT", access="g")]
def generate(self, module):
ctx = {'module': module, **self.context}
generate_file(module['module_path'], 'kernel/ExplicitEulerWithShape.templ.h', ctx)
ctx["InterfaceTestName"] = "ExplicitEulerWithShapeInterfaceCheck"
ctx["KernelInclude"] = "kernel/ExplicitEulerWithShape.h"
ctx["ExplicitInstantiation"] = "template void kernel::ExplicitEulerWithShape::operator()(const size_t p_idx1, Accessor& ac) const;"
generate_file(module['test_path'], 'tests/CheckInterface.templ.cpp', ctx,
'kernel/interfaces/ExplicitEulerWithShapeInterfaceCheck.cpp')
......@@ -5,13 +5,21 @@ from mesa_pd.utility import generate_file
class VelocityVerlet:
def __init__(self):
self.context = {'interface': []}
def __init__(self, integrate_rotation = True):
self.context = {'bIntegrateRotation': integrate_rotation, 'interface': []}
self.context['interface'].append(create_access("position", "walberla::mesa_pd::Vec3", access="gs"))
self.context['interface'].append(create_access("linearVelocity", "walberla::mesa_pd::Vec3", access="gs"))
self.context['interface'].append(create_access("invMass", "walberla::real_t", access="g"))
self.context['interface'].append(create_access("force", "walberla::mesa_pd::Vec3", access="gs"))
self.context['interface'].append(create_access("oldForce", "walberla::mesa_pd::Vec3", access="gs"))
if integrate_rotation:
self.context['interface'].append(create_access("rotation", "walberla::mesa_pd::Rot3", access="gs"))
self.context['interface'].append(create_access("angularVelocity", "walberla::mesa_pd::Vec3", access="gs"))
self.context['interface'].append(create_access("invInertiaBF", "walberla::mesa_pd::Mat3", access="g"))
self.context['interface'].append(create_access("torque", "walberla::mesa_pd::Vec3", access="gs"))
self.context['interface'].append(create_access("oldTorque", "walberla::mesa_pd::Vec3", access="gs"))
self.context['interface'].append(
create_access("flags", "walberla::mesa_pd::data::particle_flags::FlagT", access="g"))
......
# -*- coding: utf-8 -*-
from mesa_pd.accessor import create_access
from mesa_pd.utility import generate_file
class VelocityVerletWithShape:
def __init__(self):
self.context = {'interface': []}
self.context['interface'].append(create_access("position", "walberla::mesa_pd::Vec3", access="gs"))
self.context['interface'].append(create_access("linearVelocity", "walberla::mesa_pd::Vec3", access="gs"))
self.context['interface'].append(create_access("invMass", "walberla::real_t", access="g"))
self.context['interface'].append(create_access("force", "walberla::mesa_pd::Vec3", access="gs"))
self.context['interface'].append(create_access("oldForce", "walberla::mesa_pd::Vec3", access="gs"))
self.context['interface'].append(create_access("rotation", "walberla::mesa_pd::Rot3", access="gs"))
self.context['interface'].append(create_access("angularVelocity", "walberla::mesa_pd::Vec3", access="gs"))
self.context['interface'].append(create_access("invInertiaBF", "walberla::mesa_pd::Mat3", access="g"))
self.context['interface'].append(create_access("torque", "walberla::mesa_pd::Vec3", access="gs"))
self.context['interface'].append(create_access("oldTorque", "walberla::mesa_pd::Vec3", access="gs"))
self.context['interface'].append(
create_access("flags", "walberla::mesa_pd::data::particle_flags::FlagT", access="g"))
def getRequirements(self):
return self.accessor
def generate(self, module):
ctx = {'module': module, **self.context}
generate_file(module['module_path'], 'kernel/VelocityVerletWithShape.templ.h', ctx)
ctx["InterfaceTestName"] = "VelocityVerletWithShapeInterfaceCheck"
ctx["KernelInclude"] = "kernel/VelocityVerletWithShape.h"
ctx[
"ExplicitInstantiation"] = "template void kernel::VelocityVerletWithShapePreForceUpdate::operator()(const size_t p_idx1, Accessor& ac) const;\n" + \
"template void kernel::VelocityVerletWithShapePostForceUpdate::operator()(const size_t p_idx1, Accessor& ac) const;"
generate_file(module['test_path'], 'tests/CheckInterface.templ.cpp', ctx,
'kernel/interfaces/VelocityVerletWithShapeInterfaceCheck.cpp')
......@@ -2,7 +2,6 @@
from .DetectAndStoreContacts import DetectAndStoreContacts
from .DoubleCast import DoubleCast
from .ExplicitEuler import ExplicitEuler
from .ExplicitEulerWithShape import ExplicitEulerWithShape
from .ForceLJ import ForceLJ
from .HCSITSRelaxationStep import HCSITSRelaxationStep
from .HeatConduction import HeatConduction
......@@ -18,12 +17,10 @@ from .SpringDashpot import SpringDashpot
from .SpringDashpotSpring import SpringDashpotSpring
from .TemperatureIntegration import TemperatureIntegration
from .VelocityVerlet import VelocityVerlet
from .VelocityVerletWithShape import VelocityVerletWithShape
__all__ = ['DoubleCast',
'DetectAndStoreContacts',
'ExplicitEuler',
'ExplicitEulerWithShape',
'ForceLJ',
'HCSITSRelaxationStep',
'HeatConduction',
......@@ -38,5 +35,4 @@ __all__ = ['DoubleCast',
'SpringDashpot',
'SpringDashpotSpring',
'TemperatureIntegration',
'VelocityVerlet',
'VelocityVerletWithShape']
'VelocityVerlet']
......@@ -35,7 +35,6 @@ namespace kernel {
/**
* Kernel which explicitly integrates all particles in time.
* This integrator integrates velocity and position.
*
* This kernel requires the following particle accessor interface
* \code
......@@ -53,8 +52,8 @@ namespace kernel {
{%- endfor %}
* \endcode
*
* \pre All forces acting on the particles have to be set.
* \post All forces are reset to 0.
* \pre All forces and torques acting on the particles have to be set.
* \post All forces and torques are reset to 0.
* \ingroup mesa_pd_kernel
*/
class ExplicitEuler
......@@ -78,8 +77,27 @@ inline void ExplicitEuler::operator()(const size_t idx,
{
ac.setPosition (idx, ac.getInvMass(idx) * ac.getForce(idx) * dt_ * dt_ + ac.getLinearVelocity(idx) * dt_ + ac.getPosition(idx));
ac.setLinearVelocity(idx, ac.getInvMass(idx) * ac.getForce(idx) * dt_ + ac.getLinearVelocity(idx));
{%- if bIntegrateRotation %}
const Vec3 wdot = math::transformMatrixRART(ac.getRotation(idx).getMatrix(),
ac.getInvInertiaBF(idx)) * ac.getTorque(idx);
// Calculating the rotation angle
const Vec3 phi( ac.getAngularVelocity(idx) * dt_ + wdot * dt_ * dt_);
// Calculating the new orientation
auto rotation = ac.getRotation(idx);
rotation.rotate( phi );
ac.setRotation(idx, rotation);
ac.setAngularVelocity(idx, wdot * dt_ + ac.getAngularVelocity(idx));
{%- endif %}
}
ac.setForce (idx, Vec3(real_t(0), real_t(0), real_t(0)));
ac.setForce (idx, Vec3(real_t(0), real_t(0), real_t(0)));
{%- if bIntegrateRotation %}
ac.setTorque(idx, Vec3(real_t(0), real_t(0), real_t(0)));
{%- endif %}
}
} //namespace kernel
......
//======================================================================================================================
//
// This file is part of waLBerla. waLBerla is free software: you can
// redistribute it and/or modify it under the terms of the GNU General Public
// License as published by the Free Software Foundation, either version 3 of
// the License, or (at your option) any later version.
//
// waLBerla is distributed in the hope that it will be useful, but WITHOUT
// ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
// FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
// for more details.
//
// You should have received a copy of the GNU General Public License along
// with waLBerla (see COPYING.txt). If not, see <http://www.gnu.org/licenses/>.
//
//! \file ExplicitEulerWithShape.h
//! \author Sebastian Eibl <sebastian.eibl@fau.de>
//
//======================================================================================================================
//======================================================================================================================
//
// THIS FILE IS GENERATED - PLEASE CHANGE THE TEMPLATE !!!
//
//======================================================================================================================
#pragma once
#include <mesa_pd/data/DataTypes.h>
#include <mesa_pd/data/IAccessor.h>
namespace walberla {
namespace mesa_pd {
namespace kernel {
/**
* Kernel which explicitly integrates all particles in time.
* This integrator integrates velocity and position as well as angular velocity and rotation.
*
* This kernel requires the following particle accessor interface
* \code
{%- for prop in interface %}
{%- if 'g' in prop.access %}
* const {{prop.type}}& get{{prop.name | capFirst}}(const size_t p_idx) const;
{%- endif %}
{%- if 's' in prop.access %}
* void set{{prop.name | capFirst}}(const size_t p_idx, const {{prop.type}}& v);
{%- endif %}
{%- if 'r' in prop.access %}
* {{prop.type}}& get{{prop.name | capFirst}}Ref(const size_t p_idx);
{%- endif %}
*
{%- endfor %}
* \endcode
*
* \pre All forces and torques acting on the particles have to be set.
* \post All forces and torques are reset to 0.
* \ingroup mesa_pd_kernel
*/
class ExplicitEulerWithShape
{
public:
explicit ExplicitEulerWithShape(const real_t dt) : dt_(dt) {}
template <typename Accessor>
void operator()(const size_t i, Accessor& ac) const;
private:
real_t dt_ = real_t(0.0);
};
template <typename Accessor>
inline void ExplicitEulerWithShape::operator()(const size_t idx,
Accessor& ac) const
{
static_assert(std::is_base_of<data::IAccessor, Accessor>::value, "please provide a valid accessor");
if (!data::particle_flags::isSet( ac.getFlags(idx), data::particle_flags::FIXED))
{
ac.setPosition (idx, ac.getInvMass(idx) * ac.getForce(idx) * dt_ * dt_ + ac.getLinearVelocity(idx) * dt_ + ac.getPosition(idx));
ac.setLinearVelocity(idx, ac.getInvMass(idx) * ac.getForce(idx) * dt_ + ac.getLinearVelocity(idx));
const Vec3 wdot = math::transformMatrixRART(ac.getRotation(idx).getMatrix(),
ac.getInvInertiaBF(idx)) * ac.getTorque(idx);
// Calculating the rotation angle
const Vec3 phi( ac.getAngularVelocity(idx) * dt_ + wdot * dt_ * dt_);
// Calculating the new orientation
auto rotation = ac.getRotation(idx);
rotation.rotate( phi );
ac.setRotation(idx, rotation);
ac.setAngularVelocity(idx, wdot * dt_ + ac.getAngularVelocity(idx));
}
ac.setForce (idx, Vec3(real_t(0), real_t(0), real_t(0)));
ac.setTorque(idx, Vec3(real_t(0), real_t(0), real_t(0)));
}
} //namespace kernel
} //namespace mesa_pd
} //namespace walberla
......@@ -40,7 +40,6 @@ namespace kernel {
* called before the force calculation and postFroceUpdate afterwards. The
* integration is only complete when both functions are called. The integration
* is symplectic.
* \attention The force calculation has to be independent of velocity.
*
* This kernel requires the following particle accessor interface
* \code
......@@ -92,6 +91,19 @@ inline void VelocityVerletPreForceUpdate::operator()(const size_t p_idx, Accesso
ac.setPosition(p_idx, ac.getPosition(p_idx) +
ac.getLinearVelocity(p_idx) * dt_ +
real_t(0.5) * ac.getInvMass(p_idx) * ac.getOldForce(p_idx) * dt_ * dt_);
{%- if bIntegrateRotation %}
const Vec3 wdot = math::transformMatrixRART(ac.getRotation(p_idx).getMatrix(),
ac.getInvInertiaBF(p_idx)) * ac.getOldTorque(p_idx);
// Calculating the rotation angle
const Vec3 phi( ac.getAngularVelocity(p_idx) * dt_ + real_t(0.5) * wdot * dt_ * dt_);
// Calculating the new orientation
auto rotation = ac.getRotation(p_idx);
rotation.rotate( phi );
ac.setRotation(p_idx, rotation);
{%- endif %}
}
}
<