Skip to content
GitLab
Projects
Groups
Snippets
/
Help
Help
Support
Community forum
Keyboard shortcuts
?
Submit feedback
Contribute to GitLab
Sign in
Toggle navigation
Menu
Open sidebar
Dominik Mehlich
waLBerla
Commits
1f3e8215
Commit
1f3e8215
authored
Sep 16, 2019
by
Tobias Leemann
Committed by
Sebastian Eibl
Sep 16, 2019
Browse files
Integration of Hard Contact Solvers into MESA-PD
parent
1969d023
Changes
23
Expand all
Hide whitespace changes
Inline
Side-by-side
python/mesa_pd.py
View file @
1f3e8215
...
...
@@ -90,7 +90,11 @@ if __name__ == '__main__':
kernels
.
append
(
kernel
.
ExplicitEuler
()
)
kernels
.
append
(
kernel
.
ExplicitEulerWithShape
()
)
kernels
.
append
(
kernel
.
ForceLJ
()
)
kernels
.
append
(
kernel
.
HCSITSRelaxationStep
()
)
kernels
.
append
(
kernel
.
HeatConduction
()
)
kernels
.
append
(
kernel
.
InitParticlesForHCSITS
()
)
kernels
.
append
(
kernel
.
InitContactsForHCSITS
()
)
kernels
.
append
(
kernel
.
IntegrateParticlesHCSITS
()
)
kernels
.
append
(
kernel
.
InsertParticleIntoLinkedCells
()
)
kernels
.
append
(
kernel
.
LinearSpringDashpot
()
)
kernels
.
append
(
kernel
.
NonLinearSpringDashpot
()
)
...
...
python/mesa_pd/kernel/HCSITSRelaxationStep.py
0 → 100644
View file @
1f3e8215
# -*- coding: utf-8 -*-
from
mesa_pd.accessor
import
Accessor
from
..Container
import
Container
from
mesa_pd.utility
import
generateFile
class
HCSITSRelaxationStep
(
Container
):
def
__init__
(
self
):
super
().
__init__
()
self
.
addProperty
(
"maxSubIterations"
,
"size_t"
,
defValue
=
"20"
)
self
.
addProperty
(
"relaxationModel"
,
"RelaxationModel"
,
defValue
=
"InelasticFrictionlessContact"
)
self
.
addProperty
(
"deltaMax"
,
"real_t"
,
defValue
=
"0"
)
self
.
addProperty
(
"cor"
,
"real_t"
,
defValue
=
"real_t(0.2)"
)
self
.
accessor
=
Accessor
()
self
.
accessor
.
require
(
"uid"
,
"walberla::id_t"
,
access
=
"g"
)
self
.
accessor
.
require
(
"position"
,
"walberla::mesa_pd::Vec3"
,
access
=
"g"
)
self
.
accessor
.
require
(
"linearVelocity"
,
"walberla::mesa_pd::Vec3"
,
access
=
"g"
)
self
.
accessor
.
require
(
"angularVelocity"
,
"walberla::mesa_pd::Vec3"
,
access
=
"g"
)
self
.
accessor
.
require
(
"invMass"
,
"walberla::real_t"
,
access
=
"g"
)
self
.
accessor
.
require
(
"invInertia"
,
"walberla::mesa_pd::Mat3"
,
access
=
"g"
)
self
.
accessor
.
require
(
"dv"
,
"walberla::mesa_pd::Vec3"
,
access
=
"gr"
)
self
.
accessor
.
require
(
"dw"
,
"walberla::mesa_pd::Vec3"
,
access
=
"gr"
)
def
getRequirements
(
self
):
return
self
.
accessor
def
generate
(
self
,
path
):
context
=
dict
()
context
[
"properties"
]
=
self
.
properties
context
[
"interface"
]
=
self
.
accessor
.
properties
generateFile
(
path
,
'kernel/HCSITSRelaxationStep.templ.h'
,
context
)
python/mesa_pd/kernel/InitContactsForHCSITS.py
0 → 100644
View file @
1f3e8215
# -*- coding: utf-8 -*-
from
mesa_pd.accessor
import
Accessor
from
..Container
import
Container
from
mesa_pd.utility
import
generateFile
class
InitContactsForHCSITS
(
Container
):
def
__init__
(
self
):
super
().
__init__
()
self
.
addProperty
(
"erp"
,
"real_t"
,
defValue
=
"real_t(0.8)"
)
self
.
addProperty
(
"maximumPenetration"
,
"real_t"
,
defValue
=
"0"
)
self
.
paccessor
=
Accessor
()
self
.
paccessor
.
require
(
"uid"
,
"walberla::id_t"
,
access
=
"g"
)
self
.
paccessor
.
require
(
"position"
,
"walberla::mesa_pd::Vec3"
,
access
=
"g"
)
self
.
paccessor
.
require
(
"invInertia"
,
"walberla::mesa_pd::Mat3"
,
access
=
"g"
)
self
.
paccessor
.
require
(
"invMass"
,
"walberla::real_t"
,
access
=
"g"
)
def
getRequirements
(
self
):
return
self
.
paccessor
def
generate
(
self
,
path
):
context
=
dict
()
context
[
"properties"
]
=
self
.
properties
context
[
"material_parameters"
]
=
[
"friction"
]
context
[
"interface"
]
=
self
.
paccessor
.
properties
generateFile
(
path
,
'kernel/InitContactsForHCSITS.templ.h'
,
context
)
python/mesa_pd/kernel/InitParticlesForHCSITS.py
0 → 100644
View file @
1f3e8215
# -*- coding: utf-8 -*-
from
mesa_pd.accessor
import
Accessor
from
..Container
import
Container
from
mesa_pd.utility
import
generateFile
class
InitParticlesForHCSITS
(
Container
):
def
__init__
(
self
):
super
().
__init__
()
self
.
addProperty
(
"globalAcceleration"
,
"walberla::mesa_pd::Vec3"
,
defValue
=
"0"
)
self
.
accessor
=
Accessor
()
self
.
accessor
.
require
(
"uid"
,
"walberla::id_t"
,
access
=
"g"
)
self
.
accessor
.
require
(
"linearVelocity"
,
"walberla::mesa_pd::Vec3"
,
access
=
"gr"
)
self
.
accessor
.
require
(
"angularVelocity"
,
"walberla::mesa_pd::Vec3"
,
access
=
"gr"
)
self
.
accessor
.
require
(
"invMass"
,
"walberla::real_t"
,
access
=
"g"
)
self
.
accessor
.
require
(
"inertia"
,
"walberla::mesa_pd::Mat3"
,
access
=
"g"
)
self
.
accessor
.
require
(
"invInertia"
,
"walberla::mesa_pd::Mat3"
,
access
=
"g"
)
self
.
accessor
.
require
(
"dv"
,
"walberla::mesa_pd::Vec3"
,
access
=
"gr"
)
self
.
accessor
.
require
(
"dw"
,
"walberla::mesa_pd::Vec3"
,
access
=
"gr"
)
self
.
accessor
.
require
(
"torque"
,
"walberla::mesa_pd::Vec3"
,
access
=
"r"
)
self
.
accessor
.
require
(
"force"
,
"walberla::mesa_pd::Vec3"
,
access
=
"r"
)
def
getRequirements
(
self
):
return
self
.
accessor
def
generate
(
self
,
path
):
context
=
dict
()
context
[
"properties"
]
=
self
.
properties
context
[
"interface"
]
=
self
.
accessor
.
properties
generateFile
(
path
,
'kernel/InitParticlesForHCSITS.templ.h'
,
context
)
python/mesa_pd/kernel/IntegrateParticlesHCSITS.py
0 → 100644
View file @
1f3e8215
# -*- coding: utf-8 -*-
from
mesa_pd.accessor
import
Accessor
from
..Container
import
Container
from
mesa_pd.utility
import
generateFile
class
IntegrateParticlesHCSITS
(
Container
):
def
__init__
(
self
):
super
().
__init__
()
self
.
addProperty
(
"speedLimiterActive"
,
"bool"
,
defValue
=
"false"
)
self
.
addProperty
(
"speedLimitFactor"
,
"real_t"
,
defValue
=
"real_t(1.0)"
)
self
.
accessor
=
Accessor
()
self
.
accessor
.
require
(
"uid"
,
"walberla::id_t"
,
access
=
"g"
)
self
.
accessor
.
require
(
"position"
,
"walberla::mesa_pd::Vec3"
,
access
=
"gr"
)
self
.
accessor
.
require
(
"rotation"
,
"walberla::mesa_pd::Rot3"
,
access
=
"gr"
)
self
.
accessor
.
require
(
"linearVelocity"
,
"walberla::mesa_pd::Vec3"
,
access
=
"gr"
)
self
.
accessor
.
require
(
"angularVelocity"
,
"walberla::mesa_pd::Vec3"
,
access
=
"gr"
)
self
.
accessor
.
require
(
"dv"
,
"walberla::mesa_pd::Vec3"
,
access
=
"g"
)
self
.
accessor
.
require
(
"dw"
,
"walberla::mesa_pd::Vec3"
,
access
=
"g"
)
self
.
accessor
.
require
(
"flags"
,
"walberla::mesa_pd::data::particle_flags::FlagT"
,
access
=
"g"
)
def
getRequirements
(
self
):
return
self
.
accessor
def
generate
(
self
,
path
):
context
=
dict
()
context
[
"properties"
]
=
self
.
properties
context
[
"interface"
]
=
self
.
accessor
.
properties
generateFile
(
path
,
'kernel/IntegrateParticlesHCSITS.templ.h'
,
context
)
python/mesa_pd/kernel/__init__.py
View file @
1f3e8215
...
...
@@ -4,7 +4,11 @@ from .DoubleCast import DoubleCast
from
.ExplicitEuler
import
ExplicitEuler
from
.ExplicitEulerWithShape
import
ExplicitEulerWithShape
from
.ForceLJ
import
ForceLJ
from
.HCSITSRelaxationStep
import
HCSITSRelaxationStep
from
.HeatConduction
import
HeatConduction
from
.InitParticlesForHCSITS
import
InitParticlesForHCSITS
from
.InitContactsForHCSITS
import
InitContactsForHCSITS
from
.IntegrateParticlesHCSITS
import
IntegrateParticlesHCSITS
from
.InsertParticleIntoLinkedCells
import
InsertParticleIntoLinkedCells
from
.LinearSpringDashpot
import
LinearSpringDashpot
from
.NonLinearSpringDashpot
import
NonLinearSpringDashpot
...
...
@@ -19,7 +23,11 @@ __all__ = ['DoubleCast',
'ExplicitEuler'
,
'ExplicitEulerWithShape'
,
'ForceLJ'
,
'HCSITSRelaxationStep'
,
'HeatConduction'
,
'InitParticlesForHCSITS'
,
'InitContactsForHCSITS'
,
'IntegrateParticlesHCSITS'
,
'InsertParticleIntoLinkedCells'
,
'LinearSpringDashpot'
,
'NonLinearSpringDashpot'
,
...
...
python/mesa_pd/templates/data/ContactAccessor.templ.h
View file @
1f3e8215
...
...
@@ -26,7 +26,7 @@
#pragma once
#include
<mesa_pd/data/IAccessor.h>
#include
<mesa_pd/data/I
Contact
Accessor.h>
#include
<mesa_pd/data/ContactStorage.h>
#include
<core/UniqueID.h>
...
...
@@ -43,7 +43,7 @@ namespace data {
* Provides get, set and getRef for all members of the ContactStorage.
* Can be used as a basis class for a more advanced ContactAccessor.
*/
class
ContactAccessor
:
public
IAccessor
class
ContactAccessor
:
public
I
Contact
Accessor
{
public:
ContactAccessor
(
const
std
::
shared_ptr
<
data
::
ContactStorage
>&
ps
)
:
ps_
(
ps
)
{}
...
...
python/mesa_pd/templates/kernel/HCSITSRelaxationStep.templ.h
0 → 100644
View file @
1f3e8215
This diff is collapsed.
Click to expand it.
python/mesa_pd/templates/kernel/InitContactsForHCSITS.templ.h
0 → 100644
View file @
1f3e8215
//======================================================================================================================
//
// 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 InitContactsForHCSITS.h
//! \author Tobias Leemann <tobias.leemann@fau.de>
//
//======================================================================================================================
//======================================================================================================================
//
// THIS FILE IS GENERATED - PLEASE CHANGE THE TEMPLATE !!!
//
//======================================================================================================================
#pragma once
#include
<mesa_pd/common/ParticleFunctions.h>
#include
<mesa_pd/data/DataTypes.h>
#include
<mesa_pd/data/IAccessor.h>
#include
<mesa_pd/data/IContactAccessor.h>
#include
<mesa_pd/data/Flags.h>
#include
<core/logging/Logging.h>
#include
<vector>
namespace
walberla
{
namespace
mesa_pd
{
namespace
kernel
{
/**
* Init the datastructures for a contact for later use of the HCSITS-Solver.
* Call this kernel on all contacts that you want to treat with HCSITS before performing any relaxation timesteps.
* Use setErp() to set the error reduction parameter, which defines the share of the overlap that is resolved in one step
* (0 meaning after the relaxation the overlap is the same, 1 meaning the particles will have no overlap after relaxation).
* Use setFriction(a,b, cof) to define the coefficient of friction cof between materials a and b. It is assumed to be
* symmetric w.r.t. the materials.
* \ingroup mesa_pd_kernel
*/
class
InitContactsForHCSITS
{
public:
// Default constructor sets the default values
explicit
InitContactsForHCSITS
(
const
uint_t
numParticleTypes
)
:
numParticleTypes_
(
numParticleTypes
),
{
%-
for
prop
in
properties
%
}
{{
prop
.
name
}}
_
(
{{
prop
.
defValue
}}
){{
","
if
not
loop
.
last
}}
{
%-
endfor
%
}
{
{
%
for
param
in
material_parameters
%
}
{{
param
}}
_
.
resize
(
numParticleTypes
*
numParticleTypes
,
real_t
(
0
));
{
%-
endfor
%
}
}
// Getter and Setter Functions
{
%-
for
prop
in
properties
%
}
inline
const
{{
prop
.
type
}}
&
get
{{
prop
.
name
|
capFirst
}}()
const
{
return
{{
prop
.
name
}}
_
;}
inline
void
set
{{
prop
.
name
|
capFirst
}}({{
prop
.
type
}}
v
)
{
{{
prop
.
name
}}
_
=
v
;}
{
%
endfor
%
}
// Getter and Setter Functions for material parameter combinations (they are assumed symmetric).
{
%
for
param
in
material_parameters
%
}
inline
void
set
{{
param
|
capFirst
}}(
const
size_t
type1
,
const
size_t
type2
,
const
real_t
&
val
)
{
WALBERLA_ASSERT_LESS
(
type1
,
numParticleTypes_
);
WALBERLA_ASSERT_LESS
(
type2
,
numParticleTypes_
);
{{
param
}}
_
[
numParticleTypes_
*
type1
+
type2
]
=
val
;
{{
param
}}
_
[
numParticleTypes_
*
type2
+
type1
]
=
val
;
}
{
%-
endfor
%
}
{
%
for
param
in
material_parameters
%
}
inline
real_t
get
{{
param
|
capFirst
}}(
const
size_t
type1
,
const
size_t
type2
)
const
{
WALBERLA_ASSERT_LESS
(
type1
,
numParticleTypes_
);
WALBERLA_ASSERT_LESS
(
type2
,
numParticleTypes_
);
WALBERLA_ASSERT_FLOAT_EQUAL
(
{{
param
}}
_
[
numParticleTypes_
*
type1
+
type2
],
{{
param
}}
_
[
numParticleTypes_
*
type2
+
type1
],
"parameter matrix for {{param}} not symmetric!"
);
return
{{
param
}}
_
[
numParticleTypes_
*
type1
+
type2
];
}
{
%-
endfor
%
}
// List material parameters
{
%
for
param
in
material_parameters
%
}
std
::
vector
<
real_t
>
{{
param
}}
_
{};
{
%-
endfor
%
}
template
<
typename
CAccessor
,
typename
PAccessor
>
void
operator
()(
size_t
c
,
CAccessor
&
ca
,
PAccessor
&
pa
);
private:
// List properties
const
uint_t
numParticleTypes_
;
{
%
for
prop
in
properties
%
}
{{
prop
.
type
}}
{{
prop
.
name
}}
_
;
{
%-
endfor
%
}
{
%
for
param
in
parameters
%
}
std
::
vector
<
real_t
>
{{
param
}}
_
{};
{
%-
endfor
%
}
};
template
<
typename
CAccessor
,
typename
PAccessor
>
inline
void
InitContactsForHCSITS
::
operator
()(
size_t
c
,
CAccessor
&
ca
,
PAccessor
&
pa
)
{
static_assert
(
std
::
is_base_of
<
data
::
IContactAccessor
,
CAccessor
>::
value
,
"please provide a valid contact accessor"
);
static_assert
(
std
::
is_base_of
<
data
::
IAccessor
,
PAccessor
>::
value
,
"please provide a valid particle accessor"
);
size_t
bId1
=
ca
.
getId1
(
c
);
size_t
bId2
=
ca
.
getId2
(
c
);
ca
.
setR1
(
c
,
ca
.
getPosition
(
c
)
-
pa
.
getPosition
(
bId1
));
ca
.
setR2
(
c
,
ca
.
getPosition
(
c
)
-
pa
.
getPosition
(
bId2
));
// construct vector perpendicular to the normal (cross product with cardinal basis vector where the 1 component is where the other vector has its component of smallest magnitude)
if
(
std
::
fabs
(
ca
.
getNormalRef
(
c
)[
0
])
<
std
::
fabs
(
ca
.
getNormalRef
(
c
)[
1
]))
{
if
(
std
::
fabs
(
ca
.
getNormalRef
(
c
)[
0
])
<
std
::
fabs
(
ca
.
getNormalRef
(
c
)[
2
]))
ca
.
setT
(
c
,
Vec3
(
0
,
ca
.
getNormalRef
(
c
)[
2
],
-
ca
.
getNormalRef
(
c
)[
1
]));
// = n x (1, 0, 0)^T
else
ca
.
setT
(
c
,
Vec3
(
ca
.
getNormalRef
(
c
)[
1
],
-
ca
.
getNormalRef
(
c
)[
0
],
0
));
// = n x (0, 0, 1)^T
}
else
{
if
(
std
::
fabs
(
ca
.
getNormalRef
(
c
)[
1
])
<
std
::
fabs
(
ca
.
getNormalRef
(
c
)[
2
]))
ca
.
setT
(
c
,
Vec3
(
-
ca
.
getNormalRef
(
c
)[
2
],
0
,
ca
.
getNormalRef
(
c
)[
0
]));
// = n x (0, 1, 0)^T
else
ca
.
setT
(
c
,
Vec3
(
ca
.
getNormalRef
(
c
)[
1
],
-
ca
.
getNormalRef
(
c
)[
0
],
0
));
// = n x (0, 0, 1)^T
}
normalize
(
ca
.
getTRef
(
c
));
ca
.
setO
(
c
,
ca
.
getNormal
(
c
)
%
ca
.
getT
(
c
));
Mat3
contactframe
(
ca
.
getNormal
(
c
),
ca
.
getT
(
c
),
ca
.
getO
(
c
));
// If the distance is negative then penetration is present. This is an error and should be corrected.
// Correcting the whole error is not recommended since due to the linearization the errors cannot
// completely fixed anyway and the error reduction will introduce artificial restitution.
// However, if the distance is positive then it is not about error correction but the distance that
// can still be overcome without penetration and thus no error correction parameter should be applied.
if
(
ca
.
getDistance
(
c
)
<
real_t
(
0.0
))
{
setMaximumPenetration
(
std
::
max
(
getMaximumPenetration
(),
-
ca
.
getDistance
(
c
)));
ca
.
getDistanceRef
(
c
)
*=
erp_
;
}
ca
.
getMuRef
(
c
)
=
getFriction
(
pa
.
getType
(
bId1
),
pa
.
getType
(
bId2
));
Mat3
diag
=
-
(
math
::
skewSymCrossProduct
(
ca
.
getR1
(
c
),
math
::
skewSymCrossProduct
(
pa
.
getInvInertia
(
bId1
),
ca
.
getR1
(
c
)))
+
math
::
skewSymCrossProduct
(
ca
.
getR2
(
c
),
math
::
skewSymCrossProduct
(
pa
.
getInvInertia
(
bId2
),
ca
.
getR2
(
c
))));
diag
[
0
]
+=
pa
.
getInvMass
(
bId1
)
+
pa
.
getInvMass
(
bId2
);
diag
[
4
]
+=
pa
.
getInvMass
(
bId1
)
+
pa
.
getInvMass
(
bId2
);
diag
[
8
]
+=
pa
.
getInvMass
(
bId1
)
+
pa
.
getInvMass
(
bId2
);
diag
=
contactframe
.
getTranspose
()
*
diag
*
contactframe
;
// Diagonal block is known to be positive-definite and thus an inverse always exists.
ca
.
getDiag_ntoRef
(
c
)
=
diag
;
ca
.
getDiag_nto_invRef
(
c
)
=
diag
.
getInverse
();
ca
.
getDiag_n_invRef
(
c
)
=
math
::
inv
(
diag
[
0
]);
ca
.
getDiag_to_invRef
(
c
)
=
Mat2
(
diag
[
4
],
diag
[
5
],
diag
[
7
],
diag
[
8
]).
getInverse
();
}
}
}
}
python/mesa_pd/templates/kernel/InitParticlesForHCSITS.templ.h
0 → 100644
View file @
1f3e8215
//======================================================================================================================
//
// 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 InitParticlesForHCSITS.h
//! \author Tobias Leemann <tobias.leemann@fau.de>
//
//======================================================================================================================
//======================================================================================================================
//
// THIS FILE IS GENERATED - PLEASE CHANGE THE TEMPLATE !!!
//
//======================================================================================================================
#pragma once
#include
<mesa_pd/common/ParticleFunctions.h>
#include
<mesa_pd/data/DataTypes.h>
#include
<mesa_pd/data/IAccessor.h>
#include
<mesa_pd/data/Flags.h>
#include
<core/logging/Logging.h>
#include
<vector>
namespace
walberla
{
namespace
mesa_pd
{
namespace
kernel
{
/**
* Init the datastructures for the particles for later use of the HCSITS-Solver.
* Call this kernel on all particles that will be treated with HCSITS before performing any relaxation timesteps.
* Use setGlobalAcceleration() to set an accelleration action uniformly across all particles (e.g. gravity)
* \ingroup mesa_pd_kernel
*/
class
InitParticlesForHCSITS
{
public:
// Default constructor sets the default values
InitParticlesForHCSITS
()
:
{
%-
for
prop
in
properties
%
}
{{
prop
.
name
}}
_
(
{{
prop
.
defValue
}}
){{
","
if
not
loop
.
last
}}
{
%-
endfor
%
}
{}
// Getter and Setter Functions for properties
{
%-
for
prop
in
properties
%
}
inline
const
{{
prop
.
type
}}
&
get
{{
prop
.
name
|
capFirst
}}()
const
{
return
{{
prop
.
name
}}
_
;}
inline
void
set
{{
prop
.
name
|
capFirst
}}({{
prop
.
type
}}
v
)
{
{{
prop
.
name
}}
_
=
v
;}
{
%
endfor
%
}
template
<
typename
Accessor
>
void
operator
()(
size_t
j
,
Accessor
&
ac
,
real_t
dt
);
template
<
typename
Accessor
>
void
initializeVelocityCorrections
(
Accessor
&
ac
,
size_t
body
,
Vec3
&
dv
,
Vec3
&
dw
,
real_t
dt
)
const
;
private:
// List properties
{
%
for
prop
in
properties
%
}
{{
prop
.
type
}}
{{
prop
.
name
}}
_
;
{
%-
endfor
%
}
};
template
<
typename
Accessor
>
inline
void
InitParticlesForHCSITS
::
operator
()(
size_t
j
,
Accessor
&
ac
,
real_t
dt
)
{
using
namespace
data
::
particle_flags
;
static_assert
(
std
::
is_base_of
<
data
::
IAccessor
,
Accessor
>::
value
,
"please provide a valid accessor"
);
auto
particle_flags
=
ac
.
getFlagsRef
(
j
);
if
(
isSet
(
particle_flags
,
GLOBAL
)){
WALBERLA_CHECK
(
isSet
(
particle_flags
,
FIXED
),
"Global bodies need to have infinite mass as they are not communicated!"
);
initializeVelocityCorrections
(
ac
,
j
,
ac
.
getDvRef
(
j
),
ac
.
getDwRef
(
j
),
dt
);
// use applied external forces to calculate starting velocity
}
else
{
initializeVelocityCorrections
(
ac
,
j
,
ac
.
getDvRef
(
j
),
ac
.
getDwRef
(
j
),
dt
);
// use applied external forces to calculate starting velocity
if
(
!
isSet
(
particle_flags
,
FIXED
)){
// Update velocities with global acceleration and angular velocity with euler eqn
ac
.
getLinearVelocityRef
(
j
)
=
ac
.
getLinearVelocity
(
j
)
+
getGlobalAcceleration
()
*
dt
;
ac
.
getAngularVelocityRef
(
j
)
=
ac
.
getAngularVelocity
(
j
)
+
dt
*
(
ac
.
getInvInertia
(
j
)
*
(
(
ac
.
getInertia
(
j
)
*
ac
.
getAngularVelocity
(
j
)
)
%
ac
.
getAngularVelocity
(
j
)
)
);
}
}
}
//*************************************************************************************************
/*!\brief Calculates the initial velocity corrections of a given body.
* \param ac The particle accessor
* \param body The body whose velocities to time integrate
* \param dv On return the initial linear velocity correction.
* \param w On return the initial angular velocity correction.
* \param dt The time step size.
* \return void
*
* Calculates the velocity corrections effected by external forces and torques in an explicit Euler
* time integration of the velocities of the given body. For fixed objects the velocity corrections
* are set to zero. External forces and torques are reset if indicated by the settings.
*/
template
<
typename
Accessor
>
inline
void
InitParticlesForHCSITS
::
initializeVelocityCorrections
(
Accessor
&
ac
,
size_t
body
,
Vec3
&
dv
,
Vec3
&
dw
,
real_t
dt
)
const
{
dv
=
(
ac
.
getInvMass
(
body
)
*
dt
)
*
ac
.
getForce
(
body
);
dw
=
dt
*
(
ac
.
getInvInertia
(
body
)
*
ac
.
getTorque
(
body
)
);
ac
.
getForceRef
(
body
)
=
Vec3
();
ac
.
getTorqueRef
(
body
)
=
Vec3
();
}
//*************************************************************************************************
}
//namespace kernel
}
//namespace mesa_pd
}
//namespace walberla
python/mesa_pd/templates/kernel/IntegrateParticlesHCSITS.templ.h
0 → 100644
View file @
1f3e8215
//======================================================================================================================
//
// 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 IntegrateParticlesHCSITS.h
//! \author Sebastian Eibl <sebastian.eibl@fau.de>
//! \author Tobias Leemann <tobias.leemann@fau.de>
//
//======================================================================================================================
//======================================================================================================================
//
// THIS FILE IS GENERATED - PLEASE CHANGE THE TEMPLATE !!!
//
//======================================================================================================================
#pragma once
#include
<mesa_pd/common/ParticleFunctions.h>
#include
<mesa_pd/data/DataTypes.h>
#include
<mesa_pd/data/IAccessor.h>
#include
<mesa_pd/data/ContactStorage.h>
#include
<mesa_pd/data/Flags.h>
#include
<core/logging/Logging.h>
#include
<core/math/Constants.h>
#include
<vector>
namespace
walberla
{
namespace
mesa_pd
{
namespace
kernel
{
/**
* Kernel performing integration of the particles after the HCSITS iteration is done.
* Call this kernel on all particles j to integrate them by a timestep of size dt.
* The Speed Limiter limits the number of body radii, a particle can travel in one timestep.
* The speed limit factor defines a number of radii that are allowed in a timestep, e.g. a factor of 1 means, that
* a particle can only travel one time its radius in each timestep.
*
* \ingroup mesa_pd_kernel
*/
class
IntegrateParticlesHCSITS
{
public:
// Default constructor sets the default values
IntegrateParticlesHCSITS
()
:
{
%-
for
prop
in
properties
%
}
{{
prop
.
name
}}
_
(
{{
prop
.
defValue
}}
){{
","
if
not
loop
.
last
}}
{
%-
endfor
%
}
{}
// Getter and Setter Functions
{
%-
for
prop
in
properties
%
}
inline
const
{{
prop
.
type
}}
&
get
{{
prop
.
name
|
capFirst
}}()
const
{
return
{{
prop
.
name
}}
_
;}
inline
void
set
{{
prop
.
name
|
capFirst
}}({{
prop
.
type
}}
v
)
{
{{
prop
.
name
}}
_
=
v
;}
{
%
endfor
%
}
template
<
typename
PAccessor
>
void
operator
()(
size_t
j
,
PAccessor
&
ac
,
real_t
dt
);
template
<
typename
PAccessor
>
void
integratePositions
(
PAccessor
&
ac
,
size_t
body
,
Vec3
v
,
Vec3
w
,
real_t
dt
)
const
;
private:
// List properties
{
%
for
prop
in
properties
%
}
{{
prop
.
type
}}
{{
prop
.
name
}}
_
;
{
%-
endfor
%
}
};
template
<
typename
PAccessor
>
inline
void
IntegrateParticlesHCSITS
::
operator
()(
size_t
j
,
PAccessor
&
ac
,
real_t
dt
)
{
using
namespace
data
::
particle_flags
;