VelocityVerlet.templ.h 6.09 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
//======================================================================================================================
//
//  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/>.
//
16
//! \file
17
18
19
20
21
22
23
24
25
26
27
28
29
30
//! \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>
31
32
33
{%- if bIntegrateRotation %}
#include <mesa_pd/common/ParticleFunctions.h>
{%- endif %}
34
35
36
37
38
39
40
41

namespace walberla {
namespace mesa_pd {
namespace kernel {

/**
 * Velocity verlet integration for all particles.
 *
42
 * Velocity verlet integration is a two part kernel. preForceUpdate has to be
43
44
45
46
 * called before the force calculation and postFroceUpdate afterwards. The
 * integration is only complete when both functions are called. The integration
 * is symplectic.
 *
47
48
 * Wachs, A. Particle-scale computational approaches to model dry and saturated granular flows of non-Brownian, non-cohesive, and non-spherical rigid bodies. Acta Mech 230, 1919–1980 (2019). https://doi.org/10.1007/s00707-019-02389-9
 *
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
 * 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
 * \ingroup mesa_pd_kernel
 */
class VelocityVerletPreForceUpdate
{
public:
   VelocityVerletPreForceUpdate(const real_t dt) : dt_(dt) {}

   template <typename Accessor>
   void operator()(const size_t i, Accessor& ac) const;

   real_t dt_;
};

/// \see VelocityVerletPreForceUpdate
class VelocityVerletPostForceUpdate
{
public:
   VelocityVerletPostForceUpdate(const real_t dt) : dt_(dt) {}

   template <typename Accessor>
   void operator()(const size_t i, Accessor& ac) const;

   real_t dt_;
};

template <typename Accessor>
inline void VelocityVerletPreForceUpdate::operator()(const size_t p_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(p_idx), data::particle_flags::FIXED))
   {
      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_);
99
100

      {%- if bIntegrateRotation %}
101
102
103
104
105
106
107
108
      // computation done in body frame: d(omega)/ dt = J^-1 ((J*omega) x omega + T), update in world frame
      // see Wachs, 2019, doi:10.1007/s00707-019-02389-9, Eq. 27
      // note: this implementation (pre and post) is experimental as it is in principle unclear in which order
      //       angular velocities and rotations (affect again the transformations WF - BF) have to be carried out
      const auto omegaBF = transformVectorFromWFtoBF(p_idx, ac, ac.getAngularVelocity(p_idx));
      const auto torqueBF = transformVectorFromWFtoBF(p_idx, ac, ac.getOldTorque(p_idx));
      const Vec3 wdotBF = ac.getInvInertiaBF(p_idx) * ( ( ac.getInertiaBF(p_idx) * omegaBF ) % omegaBF + torqueBF );
      const Vec3 wdot = transformVectorFromBFtoWF(p_idx, ac, wdotBF);
109
110
111
112
113
114
115
116
117

      // 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 %}
118
119
120
121
122
123
124
125
126
127
128
129
   }
}

template <typename Accessor>
inline void VelocityVerletPostForceUpdate::operator()(const size_t p_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(p_idx), data::particle_flags::FIXED))
   {
      ac.setLinearVelocity(p_idx, ac.getLinearVelocity(p_idx) +
                                  real_t(0.5) * ac.getInvMass(p_idx) * (ac.getOldForce(p_idx) + ac.getForce(p_idx)) * dt_);
130
131

      {%- if bIntegrateRotation %}
132
133
134
135
      const auto omegaBF = transformVectorFromWFtoBF(p_idx, ac, ac.getAngularVelocity(p_idx));
      const auto torqueBF = transformVectorFromWFtoBF(p_idx, ac, 0.5_r * (ac.getOldTorque(p_idx) + ac.getTorque(p_idx)));
      const Vec3 wdotBF = ac.getInvInertiaBF(p_idx) * ( ( ac.getInertiaBF(p_idx) * omegaBF ) % omegaBF + torqueBF );
      const Vec3 wdot = transformVectorFromBFtoWF(p_idx, ac, wdotBF);
136
137

      ac.setAngularVelocity(p_idx, ac.getAngularVelocity(p_idx) +
138
                                   wdot * dt_ );
139
      {%- endif %}
140
   }
141

142
143
   ac.setOldForce(p_idx,       ac.getForce(p_idx));
   ac.setForce(p_idx,          Vec3(real_t(0), real_t(0), real_t(0)));
144
145
146
147
148

   {%- if bIntegrateRotation %}
   ac.setOldTorque(p_idx,      ac.getTorque(p_idx));
   ac.setTorque(p_idx,         Vec3(real_t(0), real_t(0), real_t(0)));
   {%- endif %}
149
150
151
152
153
}

} //namespace kernel
} //namespace mesa_pd
} //namespace walberla