Skip to content
GitLab
Menu
Projects
Groups
Snippets
Help
Help
Support
Community forum
Keyboard shortcuts
?
Submit feedback
Contribute to GitLab
Sign in
Toggle navigation
Menu
Open sidebar
Markus Holzer
waLBerla
Commits
1a95afc8
Commit
1a95afc8
authored
Jul 22, 2019
by
Christoph Schwarzmeier
Browse files
[API] Rename mathematical constants to make them differ from C macros
parent
57f2e7db
Changes
77
Hide whitespace changes
Inline
Side-by-side
apps/benchmarks/AdaptiveMeshRefinementFluidParticleCoupling/AMRSedimentSettling.cpp
View file @
1a95afc8
...
...
@@ -1229,7 +1229,7 @@ int main( int argc, char **argv )
const
Vector3
<
uint_t
>
domainSize
(
XBlocks
*
blockSize
,
YBlocks
*
blockSize
,
ZBlocks
*
blockSize
);
const
auto
domainVolume
=
real_t
(
domainSize
[
0
]
*
domainSize
[
1
]
*
domainSize
[
2
]);
const
real_t
sphereVolume
=
math
::
M_PI
/
real_t
(
6
)
*
diameter
*
diameter
*
diameter
;
const
real_t
sphereVolume
=
math
::
pi
/
real_t
(
6
)
*
diameter
*
diameter
*
diameter
;
const
uint_t
numberOfSediments
=
uint_c
(
std
::
ceil
(
solidVolumeFraction
*
domainVolume
/
sphereVolume
));
real_t
expectedSedimentVolumeFraction
=
(
useBox
||
useHopper
)
?
real_t
(
0.45
)
:
real_t
(
0.52
);
...
...
apps/benchmarks/AdaptiveMeshRefinementFluidParticleCoupling/WorkloadEvaluation.cpp
View file @
1a95afc8
...
...
@@ -463,8 +463,8 @@ int main( int argc, char **argv )
real_t
domainHeight
=
real_c
(
ZCells
)
-
topWallOffset
;
real_t
fluidVolume
=
real_c
(
XCells
*
YCells
)
*
domainHeight
;
real_t
solidVolume
=
solidVolumeFraction
*
fluidVolume
;
uint_t
numberOfParticles
=
uint_c
(
std
::
ceil
(
solidVolume
/
(
math
::
M_PI
/
real_t
(
6
)
*
diameter
*
diameter
*
diameter
)));
diameter
=
std
::
cbrt
(
solidVolume
/
(
real_c
(
numberOfParticles
)
*
math
::
M_PI
/
real_t
(
6
)
)
);
uint_t
numberOfParticles
=
uint_c
(
std
::
ceil
(
solidVolume
/
(
math
::
pi
/
real_t
(
6
)
*
diameter
*
diameter
*
diameter
)));
diameter
=
std
::
cbrt
(
solidVolume
/
(
real_c
(
numberOfParticles
)
*
math
::
pi
/
real_t
(
6
)
)
);
auto
densityRatio
=
real_t
(
2.5
);
...
...
@@ -821,7 +821,7 @@ int main( int argc, char **argv )
}
real_t
sphereVolume
=
diameter
*
diameter
*
diameter
*
math
::
M_PI
/
real_t
(
6
);
real_t
sphereVolume
=
diameter
*
diameter
*
diameter
*
math
::
pi
/
real_t
(
6
);
Vector3
<
real_t
>
gravitationalForce
(
real_t
(
0
),
real_t
(
0
),
-
(
densityRatio
-
real_t
(
1
))
*
gravitationalAcceleration
*
sphereVolume
);
timeloop
.
addFuncAfterTimeStep
(
pe_coupling
::
ForceOnBodiesAdder
(
blocks
,
bodyStorageID
,
gravitationalForce
),
"Gravitational force"
);
...
...
apps/benchmarks/DEM/DEM.cpp
View file @
1a95afc8
...
...
@@ -33,13 +33,13 @@ namespace dem {
real_t
calcCoefficientOfRestitution
(
const
real_t
k
,
const
real_t
gamma
,
const
real_t
meff
)
{
auto
a
=
real_t
(
0.5
)
*
gamma
/
meff
;
return
std
::
exp
(
-
a
*
math
::
M_PI
/
std
::
sqrt
(
k
/
meff
-
a
*
a
));
return
std
::
exp
(
-
a
*
math
::
pi
/
std
::
sqrt
(
k
/
meff
-
a
*
a
));
}
real_t
calcCollisionTime
(
const
real_t
k
,
const
real_t
gamma
,
const
real_t
meff
)
{
auto
a
=
real_t
(
0.5
)
*
gamma
/
meff
;
return
math
::
M_PI
/
std
::
sqrt
(
k
/
meff
-
a
*
a
);
return
math
::
pi
/
std
::
sqrt
(
k
/
meff
-
a
*
a
);
}
}
...
...
apps/benchmarks/ForcesOnSphereNearPlaneInShearFlow/ForcesOnSphereNearPlaneInShearFlow.cpp
View file @
1a95afc8
...
...
@@ -675,7 +675,7 @@ int main( int argc, char **argv )
auto
refinementTimestep
=
lbm
::
refinement
::
makeTimeStep
<
LatticeModel_T
,
BoundaryHandling_T
>
(
blocks
,
sweep
,
pdfFieldID
,
boundaryHandlingID
);
// add force evaluation and logging
real_t
normalizationFactor
=
(
zeroShearTest
)
?
real_t
(
1
)
:
(
math
::
M_PI
/
real_t
(
8
)
*
densityFluid
*
shearRate
*
shearRate
*
wallDistance
*
wallDistance
*
diameter
*
diameter
);
real_t
normalizationFactor
=
(
zeroShearTest
)
?
real_t
(
1
)
:
(
math
::
pi
/
real_t
(
8
)
*
densityFluid
*
shearRate
*
shearRate
*
wallDistance
*
wallDistance
*
diameter
*
diameter
);
std
::
string
loggingFileName
(
baseFolderLogging
+
"/LoggingForcesNearPlane"
);
loggingFileName
+=
"_lvl"
+
std
::
to_string
(
numberOfLevels
);
loggingFileName
+=
"_D"
+
std
::
to_string
(
uint_c
(
diameter
));
...
...
apps/benchmarks/MotionSingleHeavySphere/MotionSingleHeavySphere.cpp
View file @
1a95afc8
...
...
@@ -1080,7 +1080,7 @@ int main( int argc, char **argv )
WALBERLA_LOG_INFO_ON_ROOT
(
"Initial simulation has ended."
)
//evaluate the gravitational force necessary to keep the sphere at a approximately fixed position
gravity
=
forceEval
->
getForce
()
/
(
(
densityRatio
-
real_t
(
1
)
)
*
diameter
*
diameter
*
diameter
*
math
::
M_PI
/
real_t
(
6
)
);
gravity
=
forceEval
->
getForce
()
/
(
(
densityRatio
-
real_t
(
1
)
)
*
diameter
*
diameter
*
diameter
*
math
::
pi
/
real_t
(
6
)
);
GalileoSim
=
std
::
sqrt
(
(
densityRatio
-
real_t
(
1
)
)
*
gravity
*
diameter
*
diameter
*
diameter
)
/
viscosity
;
ReynoldsSim
=
uIn
*
diameter
/
viscosity
;
u_ref
=
std
::
sqrt
(
std
::
fabs
(
densityRatio
-
real_t
(
1
))
*
gravity
*
diameter
);
...
...
@@ -1234,7 +1234,7 @@ int main( int argc, char **argv )
}
// add gravity
Vector3
<
real_t
>
extForcesOnSphere
(
real_t
(
0
),
real_t
(
0
),
-
gravity
*
(
densityRatio
-
real_t
(
1
)
)
*
diameter
*
diameter
*
diameter
*
math
::
M_PI
/
real_t
(
6
));
Vector3
<
real_t
>
extForcesOnSphere
(
real_t
(
0
),
real_t
(
0
),
-
gravity
*
(
densityRatio
-
real_t
(
1
)
)
*
diameter
*
diameter
*
diameter
*
math
::
pi
/
real_t
(
6
));
timeloop
.
addFuncAfterTimeStep
(
pe_coupling
::
ForceOnBodiesAdder
(
blocks
,
bodyStorageID
,
extForcesOnSphere
),
"Add external forces (gravity and buoyancy)"
);
// evaluate the sphere properties
...
...
apps/benchmarks/PoiseuilleChannel/PoiseuilleChannel.cpp
View file @
1a95afc8
...
...
@@ -762,7 +762,7 @@ void run( const shared_ptr< Config > & config, const LatticeModel_T & latticeMod
{
setup
.
maxVelocity_L
=
(
setup
.
acceleration_L
*
setup
.
radius_L
*
setup
.
radius_L
)
/
(
real_t
(
4
)
*
setup
.
viscosity_L
);
setup
.
meanVelocity_L
=
(
setup
.
acceleration_L
*
setup
.
radius_L
*
setup
.
radius_L
)
/
(
real_t
(
8
)
*
setup
.
viscosity_L
);
setup
.
flowRate_L
=
setup
.
meanVelocity_L
*
math
::
M_PI
*
setup
.
radius_L
*
setup
.
radius_L
;
setup
.
flowRate_L
=
setup
.
meanVelocity_L
*
math
::
pi
*
setup
.
radius_L
*
setup
.
radius_L
;
}
else
{
...
...
apps/benchmarks/SchaeferTurek/SchaeferTurek.cpp
View file @
1a95afc8
...
...
@@ -744,7 +744,7 @@ public:
void
operator
()(
const
real_t
t
)
{
tConstTerm_
=
(
sinPeriod_
>
real_t
(
0
)
)
?
(
std
::
abs
(
std
::
sin
(
math
::
M_PI
*
t
/
sinPeriod_
)
)
)
:
real_t
(
1
);
tConstTerm_
=
(
sinPeriod_
>
real_t
(
0
)
)
?
(
std
::
abs
(
std
::
sin
(
math
::
pi
*
t
/
sinPeriod_
)
)
)
:
real_t
(
1
);
tConstTerm_
*=
uTerm_
*
HTerm_
;
tConstTerm_
*=
(
raisingTime_
>
real_t
(
0
)
)
?
std
::
min
(
t
/
raisingTime_
,
real_t
(
1
)
)
:
real_t
(
1
);
}
...
...
apps/showcases/BidisperseFluidizedBed/BidisperseFluidizedBedDPM.cpp
View file @
1a95afc8
...
...
@@ -371,7 +371,7 @@ uint_t createSpheresRandomly( StructuredBlockForest & forest, pe::BodyStorage &
pe
::
createSphere
(
globalBodyStorage
,
forest
.
getBlockStorage
(),
bodyStorageID
,
0
,
Vector3
<
real_t
>
(
xParticle
,
yParticle
,
zParticle
),
creationDiameter
*
real_t
(
0.5
),
material
);
currentSphereVolume
+=
math
::
M_PI
/
real_t
(
6
)
*
creationDiameter
*
creationDiameter
*
creationDiameter
;
currentSphereVolume
+=
math
::
pi
/
real_t
(
6
)
*
creationDiameter
*
creationDiameter
*
creationDiameter
;
++
numberOfSpheres
;
}
...
...
@@ -851,16 +851,16 @@ int main( int argc, char **argv ) {
const
real_t
restitutionCoeff
=
real_t
(
0.88
);
const
real_t
frictionCoeff
=
real_t
(
0.25
);
real_t
sphereVolume
=
diameterAvg
*
diameterAvg
*
diameterAvg
*
math
::
M_PI
/
real_t
(
6
);
// based on avg. diameter
real_t
sphereVolume
=
diameterAvg
*
diameterAvg
*
diameterAvg
*
math
::
pi
/
real_t
(
6
);
// based on avg. diameter
const
real_t
particleMass
=
densityRatio
*
sphereVolume
;
const
real_t
Mij
=
particleMass
*
particleMass
/
(
real_t
(
2
)
*
particleMass
);
const
real_t
lnDryResCoeff
=
std
::
log
(
restitutionCoeff
);
const
real_t
collisionTime
=
real_t
(
0.5
);
const
real_t
stiffnessCoeff
=
math
::
M_PI
*
math
::
M_PI
*
Mij
/
(
collisionTime
*
collisionTime
*
(
real_t
(
1
)
-
lnDryResCoeff
*
lnDryResCoeff
/
(
math
::
M_PI
*
math
::
M_PI
+
lnDryResCoeff
*
lnDryResCoeff
)));
const
real_t
stiffnessCoeff
=
math
::
pi
*
math
::
pi
*
Mij
/
(
collisionTime
*
collisionTime
*
(
real_t
(
1
)
-
lnDryResCoeff
*
lnDryResCoeff
/
(
math
::
pi
*
math
::
pi
+
lnDryResCoeff
*
lnDryResCoeff
)));
const
real_t
dampingCoeff
=
-
real_t
(
2
)
*
std
::
sqrt
(
Mij
*
stiffnessCoeff
)
*
(
std
::
log
(
restitutionCoeff
)
/
std
::
sqrt
(
math
::
M_PI
*
math
::
M_PI
+
(
std
::
log
(
restitutionCoeff
)
*
std
::
log
(
restitutionCoeff
))));
const
real_t
contactDuration
=
real_t
(
2
)
*
math
::
M_PI
*
Mij
/
(
std
::
sqrt
(
real_t
(
4
)
*
Mij
*
stiffnessCoeff
-
dampingCoeff
*
dampingCoeff
));
//formula from Uhlman
(
std
::
log
(
restitutionCoeff
)
/
std
::
sqrt
(
math
::
pi
*
math
::
pi
+
(
std
::
log
(
restitutionCoeff
)
*
std
::
log
(
restitutionCoeff
))));
const
real_t
contactDuration
=
real_t
(
2
)
*
math
::
pi
*
Mij
/
(
std
::
sqrt
(
real_t
(
4
)
*
Mij
*
stiffnessCoeff
-
dampingCoeff
*
dampingCoeff
));
//formula from Uhlman
WALBERLA_LOG_INFO_ON_ROOT
(
"Created particle material with:
\n
"
<<
" - coefficient of restitution = "
<<
restitutionCoeff
<<
"
\n
"
...
...
apps/tutorials/pde/01_SolvingPDE.cpp
View file @
1a95afc8
...
...
@@ -70,9 +70,9 @@ void initBC( const shared_ptr< StructuredBlockStorage > & blocks, const BlockDat
// obtain the physical coordinate of the center of the current cell
const
Vector3
<
real_t
>
p
=
blocks
->
getBlockLocalCellCenter
(
*
block
,
*
cell
);
// set the values of the Dirichlet boundary condition given by the function u(x,1) = sin(2*
M_
PI*x)*sinh(2*
M_
PI) in the source and destination field
src
->
get
(
*
cell
)
=
std
::
sin
(
real_c
(
2
)
*
math
::
M_PI
*
p
[
0
]
)
*
std
::
sinh
(
real_c
(
2
)
*
math
::
M_PI
);
dst
->
get
(
*
cell
)
=
std
::
sin
(
real_c
(
2
)
*
math
::
M_PI
*
p
[
0
]
)
*
std
::
sinh
(
real_c
(
2
)
*
math
::
M_PI
);
// set the values of the Dirichlet boundary condition given by the function u(x,1) = sin(2*PI*x)*sinh(2*PI) in the source and destination field
src
->
get
(
*
cell
)
=
std
::
sin
(
real_c
(
2
)
*
math
::
pi
*
p
[
0
]
)
*
std
::
sinh
(
real_c
(
2
)
*
math
::
pi
);
dst
->
get
(
*
cell
)
=
std
::
sin
(
real_c
(
2
)
*
math
::
pi
*
p
[
0
]
)
*
std
::
sinh
(
real_c
(
2
)
*
math
::
pi
);
}
}
}
...
...
@@ -98,9 +98,9 @@ void initRHS( const shared_ptr< StructuredBlockStorage > & blocks, const BlockDa
// obtain the physical coordinate of the center of the current cell
const
Vector3
<
real_t
>
p
=
blocks
->
getBlockLocalCellCenter
(
*
block
,
*
cell
);
// set the right-hand side, given by the function f(x,y) = 4*
M_
PI*PI*sin(2*
M_
PI*x)*sinh(2*
M_
PI*y)
f
->
get
(
*
cell
)
=
real_c
(
4
)
*
math
::
M_PI
*
math
::
M_PI
*
std
::
sin
(
real_c
(
2
)
*
math
::
M_PI
*
p
[
0
]
)
*
std
::
sinh
(
real_c
(
2
)
*
math
::
M_PI
*
p
[
1
]
);
// set the right-hand side, given by the function f(x,y) = 4*PI*PI*sin(2*PI*x)*sinh(2*PI*y)
f
->
get
(
*
cell
)
=
real_c
(
4
)
*
math
::
pi
*
math
::
pi
*
std
::
sin
(
real_c
(
2
)
*
math
::
pi
*
p
[
0
]
)
*
std
::
sinh
(
real_c
(
2
)
*
math
::
pi
*
p
[
1
]
);
}
}
}
...
...
@@ -147,7 +147,7 @@ void JacobiSweep::operator()( IBlock * const block )
dst
->
get
(
x
,
y
,
z
)
+=
(
real_c
(
1
)
/
(
dx_
*
dx_
)
)
*
src
->
get
(
x
-
1
,
y
,
z
);
dst
->
get
(
x
,
y
,
z
)
+=
(
real_c
(
1
)
/
(
dy_
*
dy_
)
)
*
src
->
get
(
x
,
y
+
1
,
z
);
dst
->
get
(
x
,
y
,
z
)
+=
(
real_c
(
1
)
/
(
dy_
*
dy_
)
)
*
src
->
get
(
x
,
y
-
1
,
z
);
dst
->
get
(
x
,
y
,
z
)
/=
(
real_c
(
2
)
/
(
dx_
*
dx_
)
+
real_c
(
2
)
/
(
dy_
*
dy_
)
+
real_c
(
4
)
*
math
::
M_PI
*
math
::
M_PI
);
dst
->
get
(
x
,
y
,
z
)
/=
(
real_c
(
2
)
/
(
dx_
*
dx_
)
+
real_c
(
2
)
/
(
dy_
*
dy_
)
+
real_c
(
4
)
*
math
::
pi
*
math
::
pi
);
)
// swap source and destination fields
...
...
@@ -282,7 +282,7 @@ int main( int argc, char ** argv )
// ...or the variant using the stencil concept
// set up the stencil weights
std
::
vector
<
real_t
>
weights
(
Stencil_T
::
Size
);
weights
[
Stencil_T
::
idx
[
stencil
::
C
]
]
=
real_c
(
2
)
/
(
dx
*
dx
)
+
real_c
(
2
)
/
(
dy
*
dy
)
+
real_c
(
4
)
*
math
::
M_PI
*
math
::
M_PI
;
weights
[
Stencil_T
::
idx
[
stencil
::
C
]
]
=
real_c
(
2
)
/
(
dx
*
dx
)
+
real_c
(
2
)
/
(
dy
*
dy
)
+
real_c
(
4
)
*
math
::
pi
*
math
::
pi
;
weights
[
Stencil_T
::
idx
[
stencil
::
E
]
]
=
real_c
(
-
1
)
/
(
dx
*
dx
);
weights
[
Stencil_T
::
idx
[
stencil
::
W
]
]
=
real_c
(
-
1
)
/
(
dx
*
dx
);
weights
[
Stencil_T
::
idx
[
stencil
::
N
]
]
=
real_c
(
-
1
)
/
(
dy
*
dy
);
...
...
apps/tutorials/pde/01_SolvingPDE.dox
View file @
1a95afc8
...
...
@@ -71,7 +71,7 @@ void JacobiSweep::operator()( IBlock * const block )
dst->get(x,y,z) += ( real_c(1) / (dx_ * dx_) ) * src->get( x-1, y , z );
dst->get(x,y,z) += ( real_c(1) / (dy_ * dy_) ) * src->get( x , y+1, z );
dst->get(x,y,z) += ( real_c(1) / (dy_ * dy_) ) * src->get( x , y-1, z );
dst->get(x,y,z) /= ( real_c(2) / (dx_ * dx_) + real_c(2)/(dy_ * dy_) + real_c(4) * math::
M_PI
* math::
M_PI
);
dst->get(x,y,z) /= ( real_c(2) / (dx_ * dx_) + real_c(2)/(dy_ * dy_) + real_c(4) * math::
pi
* math::
pi
);
)
// swap source and destination fields
...
...
@@ -90,7 +90,7 @@ typedef stencil::D2Q5 Stencil_T;
// set up the stencil weights
std::vector< real_t > weights( Stencil_T::Size );
weights[ Stencil_T::idx[ stencil::C ] ] = real_c(2) / ( dx * dx ) + real_c(2) / ( dy * dy ) + real_c(4) * math::
M_PI
* math::
M_PI
;
weights[ Stencil_T::idx[ stencil::C ] ] = real_c(2) / ( dx * dx ) + real_c(2) / ( dy * dy ) + real_c(4) * math::
pi
* math::
pi
;
weights[ Stencil_T::idx[ stencil::E ] ] = real_c(-1) / ( dx * dx );
weights[ Stencil_T::idx[ stencil::W ] ] = real_c(-1) / ( dx * dx );
weights[ Stencil_T::idx[ stencil::N ] ] = real_c(-1) / ( dy * dy );
...
...
@@ -185,9 +185,9 @@ void initRHS( const shared_ptr< StructuredBlockStorage > & blocks, const BlockDa
// obtain the physical coordinate of the center of the current cell
const Vector3< real_t > p = blocks->getBlockLocalCellCenter( *block, *cell );
// set the right-hand side, given by the function f(x,y) = 4*
M_
PI*PI*sin(2*
M_
PI*x)*sinh(2*
M_
PI*y)
f->get( *cell ) = real_c(4) * math::
M_PI
* math::
M_PI
* std::sin( real_c(2) * math::
M_PI
* p[0] ) *
std::sinh( real_c(2) * math::
M_PI
* p[1] );
// set the right-hand side, given by the function f(x,y) = 4*PI*PI*sin(2*PI*x)*sinh(2*PI*y)
f->get( *cell ) = real_c(4) * math::
pi
* math::
pi
* std::sin( real_c(2) * math::
pi
* p[0] ) *
std::sinh( real_c(2) * math::
pi
* p[1] );
}
}
}
...
...
@@ -236,10 +236,10 @@ void initBC( const shared_ptr< StructuredBlockStorage > & blocks, const BlockDat
// obtain the physical coordinate of the center of the current cell
const Vector3< real_t > p = blocks->getBlockLocalCellCenter( *block, *cell );
// set the values of the Dirichlet boundary condition given by the function u(x,1) = sin(2*
M_
PI*x)*sinh(2*
M_
PI)
// set the values of the Dirichlet boundary condition given by the function u(x,1) = sin(2*PI*x)*sinh(2*PI)
// in the source and destination field
src->get( *cell ) = std::sin( real_c(2) * math::
M_PI
* p[0] ) * std::sinh( real_c(2) * math::
M_PI
);
dst->get( *cell ) = std::sin( real_c(2) * math::
M_PI
* p[0] ) * std::sinh( real_c(2) * math::
M_PI
);
src->get( *cell ) = std::sin( real_c(2) * math::
pi
* p[0] ) * std::sinh( real_c(2) * math::
pi
);
dst->get( *cell ) = std::sin( real_c(2) * math::
pi
* p[0] ) * std::sinh( real_c(2) * math::
pi
);
}
}
}
...
...
apps/tutorials/pde/02_HeatEquation.cpp
View file @
1a95afc8
...
...
@@ -66,8 +66,8 @@ void initU( const shared_ptr< StructuredBlockStorage > & blocks, const BlockData
// obtain the physical coordinate of the center of the current cell
const
Vector3
<
real_t
>
p
=
blocks
->
getBlockLocalCellCenter
(
*
block
,
*
cell
);
// set the initial condition, given by the function u(x,y,0) = sin(
M_
PI*x)*sin(
M_
PI*y)
u
->
get
(
*
cell
)
=
std
::
sin
(
math
::
M_PI
*
p
[
0
]
)
*
std
::
sin
(
math
::
M_PI
*
p
[
1
]
);
// set the initial condition, given by the function u(x,y,0) = sin(PI*x)*sin(PI*y)
u
->
get
(
*
cell
)
=
std
::
sin
(
math
::
pi
*
p
[
0
]
)
*
std
::
sin
(
math
::
pi
*
p
[
1
]
);
}
}
}
...
...
apps/tutorials/pde/03_HeatEquation_Extensions.cpp
View file @
1a95afc8
...
...
@@ -67,8 +67,8 @@ void initU( const shared_ptr< StructuredBlockStorage > & blocks, const BlockData
// obtain the physical coordinate of the center of the current cell
const
Vector3
<
real_t
>
p
=
blocks
->
getBlockLocalCellCenter
(
*
block
,
*
cell
);
// set the initial condition, given by the function u(x,y,0) = sin(
M_
PI*x)*sin(
M_
PI*y)
u
->
get
(
*
cell
)
=
std
::
sin
(
math
::
M_PI
*
p
[
0
]
)
*
std
::
sin
(
math
::
M_PI
*
p
[
1
]
);
// set the initial condition, given by the function u(x,y,0) = sin(PI*x)*sin(PI*y)
u
->
get
(
*
cell
)
=
std
::
sin
(
math
::
pi
*
p
[
0
]
)
*
std
::
sin
(
math
::
pi
*
p
[
1
]
);
}
}
}
...
...
python/mesa_pd/templates/kernel/SpringDashpot.templ.h
View file @
1a95afc8
...
...
@@ -90,7 +90,7 @@ public:
const
real_t
meff
)
{
auto
a
=
real_t
(
0.5
)
*
getDampingN
(
type1
,
type2
)
/
meff
;
return
std
::
exp
(
-
a
*
math
::
M_PI
/
std
::
sqrt
(
getStiffness
(
type1
,
type2
)
/
meff
-
a
*
a
));
return
std
::
exp
(
-
a
*
math
::
pi
/
std
::
sqrt
(
getStiffness
(
type1
,
type2
)
/
meff
-
a
*
a
));
}
inline
...
...
@@ -99,7 +99,7 @@ public:
const
real_t
meff
)
{
auto
a
=
real_t
(
0.5
)
*
getDampingN
(
type1
,
type2
)
/
meff
;
return
math
::
M_PI
/
std
::
sqrt
(
getStiffness
(
type1
,
type2
)
/
meff
-
a
*
a
);
return
math
::
pi
/
std
::
sqrt
(
getStiffness
(
type1
,
type2
)
/
meff
-
a
*
a
);
}
inline
...
...
@@ -110,8 +110,8 @@ public:
const
real_t
meff
)
{
const
real_t
lnDryResCoeff
=
std
::
log
(
cor
);
setStiffness
(
type1
,
type2
,
math
::
M_PI
*
math
::
M_PI
*
meff
/
(
collisionTime
*
collisionTime
*
(
real_t
(
1
)
-
lnDryResCoeff
*
lnDryResCoeff
/
(
math
::
M_PI
*
math
::
M_PI
+
lnDryResCoeff
*
lnDryResCoeff
))
));
setDampingN
(
type1
,
type2
,
-
real_t
(
2
)
*
std
::
sqrt
(
meff
*
getStiffness
(
type1
,
type2
)
)
*
(
lnDryResCoeff
/
std
::
sqrt
(
math
::
M_PI
*
math
::
M_PI
+
(
lnDryResCoeff
*
lnDryResCoeff
)
)
));
setStiffness
(
type1
,
type2
,
math
::
pi
*
math
::
pi
*
meff
/
(
collisionTime
*
collisionTime
*
(
real_t
(
1
)
-
lnDryResCoeff
*
lnDryResCoeff
/
(
math
::
pi
*
math
::
pi
+
lnDryResCoeff
*
lnDryResCoeff
))
));
setDampingN
(
type1
,
type2
,
-
real_t
(
2
)
*
std
::
sqrt
(
meff
*
getStiffness
(
type1
,
type2
)
)
*
(
lnDryResCoeff
/
std
::
sqrt
(
math
::
pi
*
math
::
pi
+
(
lnDryResCoeff
*
lnDryResCoeff
)
)
));
}
private:
uint_t
numParticleTypes_
;
...
...
src/core/math/Constants.h
View file @
1a95afc8
...
...
@@ -16,57 +16,43 @@
//! \file Constants.h
//! \author Klaus Iglberger
//! \author Sebastian Eibl <sebastian.eibl@fau.de>
//! \author Christoph Schwarzmeier <christoph.schwarzmeier@fau.de>
//! \brief Header file for mathematical constants
//
//======================================================================================================================
#pragma once
//*************************************************************************************************
// Includes
//*************************************************************************************************
#include <core/DataTypes.h>
#include <cmath>
#include <core/DataTypes.h>
// Disable false warnings in GCC 5
#if (
defined __GNUC__
) && (
__GNUC__ == 5
) && (
__GNUC_MINOR__ == 1
)
# pragma GCC diagnostic push
# pragma GCC diagnostic ignored "-Wunused-variable"
#if (defined __GNUC__) && (__GNUC__ == 5) && (__GNUC_MINOR__ == 1)
#
pragma GCC diagnostic push
#
pragma GCC diagnostic ignored "-Wunused-variable"
#endif
namespace
walberla
{
namespace
math
{
# undef M_E
# undef M_LOG2E
# undef M_LOG10E
# undef M_LN2
# undef M_LN10
# undef M_PI
# undef M_PI_2
# undef M_PI_4
# undef M_1_PI
# undef M_2_PI
# undef M_2_SQRTPI
# undef M_SQRT2
# undef M_SQRT1_2
constexpr
real_t
M_E
=
real_t
(
2.7182818284590452354
);
/* e */
constexpr
real_t
M_LOG2E
=
real_t
(
1.4426950408889634074
);
/* log_2 e */
constexpr
real_t
M_LOG10E
=
real_t
(
0.43429448190325182765
);
/* log_10 e */
constexpr
real_t
M_LN2
=
real_t
(
0.69314718055994530942
);
/* log_e 2 */
constexpr
real_t
M_LN10
=
real_t
(
2.30258509299404568402
);
/* log_e 10 */
constexpr
real_t
M_PI
=
real_t
(
3.14159265358979323846
);
/* pi */
constexpr
real_t
M_PI_2
=
real_t
(
1.57079632679489661923
);
/* pi/2 */
constexpr
real_t
M_PI_4
=
real_t
(
0.78539816339744830962
);
/* pi/4 */
constexpr
real_t
M_1_PI
=
real_t
(
0.31830988618379067154
);
/* 1/pi */
constexpr
real_t
M_2_PI
=
real_t
(
0.63661977236758134308
);
/* 2/pi */
constexpr
real_t
M_2_SQRTPI
=
real_t
(
1.12837916709551257390
);
/* 2/sqrt(pi) */
constexpr
real_t
M_SQRT2
=
real_t
(
1.41421356237309504880
);
/* sqrt(2) */
constexpr
real_t
M_SQRT1_2
=
real_t
(
0.70710678118654752440
);
/* 1/sqrt(2) */
namespace
walberla
{
namespace
math
{
constexpr
real_t
e
=
real_t
(
2.7182818284590452354
);
/* e */
constexpr
real_t
log2_e
=
real_t
(
1.4426950408889634074
);
/* log_2 e */
constexpr
real_t
log10_e
=
real_t
(
0.43429448190325182765
);
/* log_10 e */
constexpr
real_t
ln_two
=
real_t
(
0.69314718055994530942
);
/* log_e 2 */
constexpr
real_t
ln_ten
=
real_t
(
2.30258509299404568402
);
/* log_e 10 */
constexpr
real_t
pi
=
real_t
(
3.14159265358979323846
);
/* pi */
constexpr
real_t
half_pi
=
real_t
(
1.57079632679489661923
);
/* pi/2 */
constexpr
real_t
fourth_pi
=
real_t
(
0.78539816339744830962
);
/* pi/4 */
constexpr
real_t
one_div_pi
=
real_t
(
0.31830988618379067154
);
/* 1/pi */
constexpr
real_t
two_div_pi
=
real_t
(
0.63661977236758134308
);
/* 2/pi */
constexpr
real_t
two_div_root_pi
=
real_t
(
1.12837916709551257390
);
/* 2/sqrt(pi) */
constexpr
real_t
root_two
=
real_t
(
1.41421356237309504880
);
/* sqrt(2) */
constexpr
real_t
one_div_root_two
=
real_t
(
0.70710678118654752440
);
/* 1/sqrt(2) */
}
// namespace math
}
}
// namespace walberla
\ No newline at end of file
src/core/math/GenericAABB.h
View file @
1a95afc8
...
...
@@ -111,8 +111,8 @@ public:
inline
bool
containsClosedInterval
(
const
vector_type
&
point
)
const
;
inline
bool
containsClosedInterval
(
const
vector_type
&
point
,
const
value_type
dx
)
const
;
inline
GenericAABB
getExtended
(
const
value_type
e
)
const
;
inline
GenericAABB
getExtended
(
const
vector_type
&
e
)
const
;
inline
GenericAABB
getExtended
(
const
value_type
e
ps
)
const
;
inline
GenericAABB
getExtended
(
const
vector_type
&
e
ps
)
const
;
inline
GenericAABB
getTranslated
(
const
vector_type
&
translation
)
const
;
...
...
@@ -162,8 +162,8 @@ public:
inline
void
setAxisBounds
(
const
uint_t
index
,
const
value_type
value1
,
const
value_type
value2
);
inline
void
extend
(
const
value_type
e
);
inline
void
extend
(
const
vector_type
&
e
);
inline
void
extend
(
const
value_type
e
ps
);
inline
void
extend
(
const
vector_type
&
e
ps
);
inline
void
setCenter
(
const
vector_type
&
center
);
inline
void
translate
(
const
vector_type
&
translation
);
...
...
src/core/math/GenericAABB.impl.h
View file @
1a95afc8
...
...
@@ -674,19 +674,19 @@ bool GenericAABB< T >::containsClosedInterval( const vector_type & point, const
/**
* \brief Creates a new GenericAABB by extending this one
*
* \param e epsilon by which the bounding box is extended in each direction
* \param e
ps
epsilon by which the bounding box is extended in each direction
*
* \returns The extended GenericAABB
*/
template
<
typename
T
>
GenericAABB
<
T
>
GenericAABB
<
T
>::
getExtended
(
const
value_type
e
)
const
GenericAABB
<
T
>
GenericAABB
<
T
>::
getExtended
(
const
value_type
e
ps
)
const
{
vector_type
newMinCorner
(
minCorner_
[
0
]
-
e
,
minCorner_
[
1
]
-
e
,
minCorner_
[
2
]
-
e
);
vector_type
newMinCorner
(
minCorner_
[
0
]
-
e
ps
,
minCorner_
[
1
]
-
e
ps
,
minCorner_
[
2
]
-
e
ps
);
return
createFromMinMaxCorner
(
newMinCorner
[
0
],
newMinCorner
[
1
],
newMinCorner
[
2
],
std
::
max
(
newMinCorner
[
0
],
maxCorner_
[
0
]
+
e
),
std
::
max
(
newMinCorner
[
1
],
maxCorner_
[
1
]
+
e
),
std
::
max
(
newMinCorner
[
2
],
maxCorner_
[
2
]
+
e
)
);
std
::
max
(
newMinCorner
[
0
],
maxCorner_
[
0
]
+
e
ps
),
std
::
max
(
newMinCorner
[
1
],
maxCorner_
[
1
]
+
e
ps
),
std
::
max
(
newMinCorner
[
2
],
maxCorner_
[
2
]
+
e
ps
)
);
}
...
...
@@ -694,19 +694,19 @@ GenericAABB< T > GenericAABB< T >::getExtended( const value_type e ) const
/**
* \brief Creates a new GenericAABB by extending this one
*
* \param e epsilon vector by which the bounding box is extended. The box is extended in each direction by
* \param e
ps
epsilon vector by which the bounding box is extended. The box is extended in each direction by
* the corresponding vector component.
*
* \returns The extended GenericAABB
*/
template
<
typename
T
>
GenericAABB
<
T
>
GenericAABB
<
T
>::
getExtended
(
const
vector_type
&
e
)
const
GenericAABB
<
T
>
GenericAABB
<
T
>::
getExtended
(
const
vector_type
&
e
ps
)
const
{
vector_type
newMinCorner
(
minCorner_
-
e
);
vector_type
newMinCorner
(
minCorner_
-
e
ps
);
return
createFromMinMaxCorner
(
newMinCorner
[
0
],
newMinCorner
[
1
],
newMinCorner
[
2
],
std
::
max
(
newMinCorner
[
0
],
maxCorner_
[
0
]
+
e
[
0
]
),
std
::
max
(
newMinCorner
[
1
],
maxCorner_
[
1
]
+
e
[
0
]
),
std
::
max
(
newMinCorner
[
2
],
maxCorner_
[
2
]
+
e
[
0
]
)
);
std
::
max
(
newMinCorner
[
0
],
maxCorner_
[
0
]
+
e
ps
[
0
]
),
std
::
max
(
newMinCorner
[
1
],
maxCorner_
[
1
]
+
e
ps
[
0
]
),
std
::
max
(
newMinCorner
[
2
],
maxCorner_
[
2
]
+
e
ps
[
0
]
)
);
}
...
...
@@ -1519,18 +1519,18 @@ void GenericAABB< T >::setAxisBounds( const uint_t index, const value_type value
/**
* \brief Extends this GenericAABB
*
* \param e epsilon by which the bounding box is extended in each direction
* \param e
ps
epsilon by which the bounding box is extended in each direction
*/
template
<
typename
T
>
void
GenericAABB
<
T
>::
extend
(
const
value_type
e
)
void
GenericAABB
<
T
>::
extend
(
const
value_type
e
ps
)
{
minCorner_
[
0
]
-=
e
;
minCorner_
[
1
]
-=
e
;
minCorner_
[
2
]
-=
e
;
minCorner_
[
0
]
-=
e
ps
;
minCorner_
[
1
]
-=
e
ps
;
minCorner_
[
2
]
-=
e
ps
;
maxCorner_
[
0
]
+=
e
;
maxCorner_
[
1
]
+=
e
;
maxCorner_
[
2
]
+=
e
;
maxCorner_
[
0
]
+=
e
ps
;
maxCorner_
[
1
]
+=
e
ps
;
maxCorner_
[
2
]
+=
e
ps
;
for
(
uint_t
i
=
0
;
i
<
3
;
++
i
)
if
(
minCorner_
[
i
]
>
maxCorner_
[
i
]
)
...
...
@@ -1544,14 +1544,14 @@ void GenericAABB< T >::extend( const value_type e )
/**
* \brief Extends this GenericAABB
*
* \param e epsilon vector by which the bounding box is extended. The box is extended in each direction by
* \param e
ps
epsilon vector by which the bounding box is extended. The box is extended in each direction by
* the corresponding vector component.
*/
template
<
typename
T
>
void
GenericAABB
<
T
>::
extend
(
const
vector_type
&
e
)
void
GenericAABB
<
T
>::
extend
(
const
vector_type
&
e
ps
)
{
minCorner_
-=
e
;
maxCorner_
+=
e
;
minCorner_
-=
e
ps
;
maxCorner_
+=
e
ps
;
for
(
uint_t
i
=
0
;
i
<
3
;
++
i
)
if
(
minCorner_
[
i
]
>
maxCorner_
[
i
]
)
...
...
src/core/math/equation_system/EquationParser.cpp
View file @
1a95afc8
...
...
@@ -189,12 +189,12 @@ NodePtr EquationParser::parseFunction( const std::string& str, size_t& index ) c
{
case
OP_FUNC_EXP
:
funcPtr
=
std
::
make_shared
<
Node
>
(
OP_PROD
);
funcPtr
->
left
()
=
std
::
make_shared
<
Node
>
(
M_E
);
funcPtr
->
left
()
=
std
::
make_shared
<
Node
>
(
math
::
e
);
funcPtr
->
right
()
=
nodePtr
;
return
funcPtr
;
case
OP_FUNC_LN
:
funcPtr
=
std
::
make_shared
<
Node
>
(
OP_LOG
);
funcPtr
->
right
()
=
std
::
make_shared
<
Node
>
(
M_E
);
funcPtr
->
right
()
=
std
::
make_shared
<
Node
>
(
math
::
e
);
funcPtr
->
left
()
=
nodePtr
;
return
funcPtr
;
case
OP_FUNC_SQRT
:
...
...
src/mesa_pd/data/shape/Sphere.h
View file @
1a95afc8
...
...
@@ -37,7 +37,7 @@ public:
void
updateMassAndInertia
(
const
real_t
density
)
override
;
real_t
getVolume
()
const
override
{
return
(
real_t
(
4
)
/
real_t
(
3
))
*
math
::
M_PI
*
getRadius
()
*
getRadius
()
*
getRadius
();
}
real_t
getVolume
()
const
override
{
return
(
real_t
(
4
)
/
real_t
(
3
))
*
math
::
pi
*
getRadius
()
*
getRadius
()
*
getRadius
();
}
static
const
int
SHAPE_TYPE
=
1
;
///< Unique shape type identifier for spheres.\ingroup mesa_pd_shape
...
...
@@ -48,7 +48,7 @@ private:
inline
void
Sphere
::
updateMassAndInertia
(
const
real_t
density
)
{
const
real_t
m
=
(
real_c
(
4.0
)
/
real_c
(
3.0
)
*
math
::
M_PI
)
*
getRadius
()
*
getRadius
()
*
getRadius
()
*
density
;
const
real_t
m
=
(
real_c
(
4.0
)
/
real_c
(
3.0
)
*
math
::
pi
)
*
getRadius
()
*
getRadius
()
*
getRadius
()
*
density
;
const
Mat3
I
=
Mat3
::
makeDiagonalMatrix
(
real_c
(
0.4
)
*
m
*
getRadius
()
*
getRadius
()
);
mass_
=
m
;
...
...
src/mesa_pd/kernel/SpringDashpot.h
View file @
1a95afc8
...
...
@@ -100,7 +100,7 @@ public:
const
real_t
meff
)
{
auto
a
=
real_t
(
0.5
)
*
getDampingN
(
type1
,
type2
)
/
meff
;
return
std
::
exp
(
-
a
*
math
::
M_PI
/
std
::
sqrt
(
getStiffness
(
type1
,
type2
)
/
meff
-
a
*
a
));
return
std
::
exp
(
-
a
*
math
::
pi
/
std
::
sqrt
(
getStiffness
(
type1
,
type2
)
/
meff
-
a
*
a
));
}
inline
...
...
@@ -109,7 +109,7 @@ public:
const
real_t
meff
)
{
auto
a
=
real_t
(
0.5
)
*
getDampingN
(
type1
,
type2
)
/
meff
;
return
math
::
M_PI
/
std
::
sqrt
(
getStiffness
(
type1
,
type2
)
/
meff
-
a
*
a
);
return
math
::
pi
/
std
::
sqrt
(
getStiffness
(
type1
,
type2
)
/
meff
-
a
*
a
);
}
inline
...
...
@@ -120,8 +120,8 @@ public:
const
real_t
meff
)
{
const
real_t
lnDryResCoeff
=
std
::
log
(
cor
);
setStiffness
(
type1
,
type2
,
math
::
M_PI
*
math
::
M_PI
*
meff
/
(
collisionTime
*
collisionTime
*
(
real_t
(
1
)
-
lnDryResCoeff
*
lnDryResCoeff
/
(
math
::
M_PI
*
math
::
M_PI
+
lnDryResCoeff
*
lnDryResCoeff
))
));
setDampingN
(
type1
,
type2
,
-
real_t
(
2
)
*
std
::
sqrt
(
meff
*
getStiffness
(
type1
,
type2
)
)
*
(
lnDryResCoeff
/
std
::
sqrt
(
math
::
M_PI
*
math
::
M_PI
+
(
lnDryResCoeff
*
lnDryResCoeff
)
)
));
setStiffness
(
type1
,
type2
,
math
::
pi
*
math
::
pi
*
meff
/
(
collisionTime
*
collisionTime
*
(
real_t
(
1
)
-
lnDryResCoeff
*
lnDryResCoeff
/
(
math
::
pi
*
math
::
pi
+
lnDryResCoeff
*
lnDryResCoeff
))
));
setDampingN
(
type1
,
type2
,
-
real_t
(
2
)
*
std
::
sqrt
(
meff
*
getStiffness
(
type1
,
type2
)
)
*
(
lnDryResCoeff
/
std
::
sqrt
(
math
::
pi
*
math
::
pi
+
(
lnDryResCoeff
*
lnDryResCoeff
)
)
));
}
private:
uint_t
numParticleTypes_
;
...
...
src/pe/collision/EPA.cpp
View file @
1a95afc8
...
...
@@ -536,7 +536,7 @@ inline void EPA::createInitialSimplex( size_t numPoints, GeomPrimitive &geom1, G
}
Vec3
direction1
=
(
d
%
axis
).
getNormalized
();
Quat
q
(
d
,
(
real_t
(
2.0
)
/
real_t
(
3.0
))
*
real_t
(
walberla
::
math
::
M_PI
));
Quat
q
(
d
,
(
real_t
(
2.0
)
/
real_t
(
3.0
))
*
real_t
(
walberla
::
math
::
pi
));
Mat3
rot
=
q
.
toRotationMatrix
();
Vec3
direction2
=
(
rot
*
direction1
).
getNormalized
();
Vec3
direction3
=
(
rot
*
direction2
).
getNormalized
();
...
...
Prev
1
2
3
4
Next
Write
Preview
Markdown
is supported
0%
Try again
or
attach a new file
.
Attach a file
Cancel
You are about to add
0
people
to the discussion. Proceed with caution.
Finish editing this message first!
Cancel
Please
register
or
sign in
to comment