Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F68302120
test_solid_mechanics_model_dynamics.cc
No One
Temporary
Actions
Download File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Subscribers
None
File Metadata
Details
File Info
Storage
Attached
Created
Wed, Jun 26, 19:40
Size
8 KB
Mime Type
text/x-c++
Expires
Fri, Jun 28, 19:40 (2 d)
Engine
blob
Format
Raw Data
Handle
18378671
Attached To
rAKA akantu
test_solid_mechanics_model_dynamics.cc
View Options
/**
* @file test_solid_mechanics_model_cube3d.cc
*
* @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
*
* @date creation: Wed Aug 04 2010
* @date last modification: Thu Aug 06 2015
*
* @brief test of the class SolidMechanicsModel on the 3d cube
*
* @section LICENSE
*
* Copyright (©) 2010-2012, 2014, 2015 EPFL (Ecole Polytechnique Fédérale de
* Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des
* Solides)
*
* Akantu is free software: you can redistribute it and/or modify it under the
* terms of the GNU Lesser General Public License as published by the Free
* Software Foundation, either version 3 of the License, or (at your option) any
* later version.
*
* Akantu 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 Lesser General Public License for more
* details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with Akantu. If not, see <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#include "boundary_condition_functor.hh"
#include "test_solid_mechanics_model_fixture.hh"
/* -------------------------------------------------------------------------- */
using
namespace
akantu
;
namespace
{
constexpr
Real
A
=
1e-1
;
constexpr
Real
E
=
1.
;
constexpr
Real
poisson
=
.3
;
constexpr
Real
lambda
=
E
*
poisson
/
(
1.
+
poisson
)
/
(
1.
-
2.
*
poisson
);
constexpr
Real
mu
=
E
/
2
/
(
1.
+
poisson
);
constexpr
Real
rho
=
1.
;
const
Real
cp
=
std
::
sqrt
((
lambda
+
2
*
mu
)
/
rho
);
const
Real
cs
=
std
::
sqrt
(
mu
/
rho
);
const
Real
c
=
std
::
sqrt
(
E
/
rho
);
const
Vector
<
Real
>
k
=
{
.5
,
0.
,
0.
};
const
Vector
<
Real
>
psi1
=
{
0.
,
0.
,
1.
};
const
Vector
<
Real
>
psi2
=
{
0.
,
1.
,
0.
};
const
Real
knorm
=
k
.
norm
();
template
<
UInt
spatial_dimension
>
auto
solution_disp
=
1
;
template
<
UInt
spatial_dimension
>
auto
solution_vel
=
1
;
template
<>
auto
solution_disp
<
1
>
=
[](
Vector
<
Real
>
&
disp
,
const
Vector
<
Real
>
&
coord
,
Real
current_time
)
{
const
auto
&
x
=
coord
(
_x
);
const
Real
omega
=
c
*
k
[
0
];
disp
(
_x
)
=
A
*
cos
(
k
[
0
]
*
x
-
omega
*
current_time
);
};
template
<>
auto
solution_vel
<
1
>
=
[](
Vector
<
Real
>
&
vel
,
const
Vector
<
Real
>
&
coord
,
Real
current_time
)
{
const
auto
&
x
=
coord
(
_x
);
const
Real
omega
=
c
*
k
[
0
];
vel
(
_x
)
=
omega
*
A
*
sin
(
k
[
0
]
*
x
-
omega
*
current_time
);
};
template
<>
auto
solution_disp
<
2
>
=
[](
Vector
<
Real
>
&
disp
,
const
Vector
<
Real
>
&
coord
,
Real
current_time
)
{
const
auto
&
X
=
coord
;
Vector
<
Real
>
kshear
=
{
k
[
1
],
k
[
0
]};
Vector
<
Real
>
kpush
=
{
k
[
0
],
k
[
1
]};
const
Real
omega_p
=
knorm
*
cp
;
const
Real
omega_s
=
knorm
*
cs
;
Real
phase_p
=
X
.
dot
(
kpush
)
-
omega_p
*
current_time
;
Real
phase_s
=
X
.
dot
(
kpush
)
-
omega_s
*
current_time
;
disp
=
A
*
(
kpush
*
cos
(
phase_p
)
+
kshear
*
cos
(
phase_s
));
};
template
<>
auto
solution_vel
<
2
>
=
[](
Vector
<
Real
>
&
vel
,
const
Vector
<
Real
>
&
coord
,
Real
current_time
)
{
const
auto
&
X
=
coord
;
Vector
<
Real
>
kshear
=
{
k
[
1
],
k
[
0
]};
Vector
<
Real
>
kpush
=
{
k
[
0
],
k
[
1
]};
const
Real
omega_p
=
knorm
*
cp
;
const
Real
omega_s
=
knorm
*
cs
;
Real
phase_p
=
X
.
dot
(
kpush
)
-
omega_p
*
current_time
;
Real
phase_s
=
X
.
dot
(
kpush
)
-
omega_s
*
current_time
;
vel
=
A
*
(
kpush
*
omega_p
*
cos
(
phase_p
)
+
kshear
*
omega_s
*
cos
(
phase_s
));
};
template
<>
auto
solution_disp
<
3
>
=
[](
Vector
<
Real
>
&
disp
,
const
Vector
<
Real
>
&
coord
,
Real
current_time
)
{
const
auto
&
X
=
coord
;
Vector
<
Real
>
kpush
=
k
;
Vector
<
Real
>
kshear1
(
3
);
Vector
<
Real
>
kshear2
(
3
);
kshear1
.
crossProduct
(
k
,
psi1
);
kshear2
.
crossProduct
(
k
,
psi2
);
const
Real
omega_p
=
knorm
*
cp
;
const
Real
omega_s
=
knorm
*
cs
;
Real
phase_p
=
X
.
dot
(
kpush
)
-
omega_p
*
current_time
;
Real
phase_s
=
X
.
dot
(
kpush
)
-
omega_s
*
current_time
;
disp
=
A
*
(
kpush
*
cos
(
phase_p
)
+
kshear1
*
cos
(
phase_s
)
+
kshear2
*
cos
(
phase_s
));
};
template
<>
auto
solution_vel
<
3
>
=
[](
Vector
<
Real
>
&
vel
,
const
Vector
<
Real
>
&
coord
,
Real
current_time
)
{
const
auto
&
X
=
coord
;
Vector
<
Real
>
kpush
=
k
;
Vector
<
Real
>
kshear1
(
3
);
Vector
<
Real
>
kshear2
(
3
);
kshear1
.
crossProduct
(
k
,
psi1
);
kshear2
.
crossProduct
(
k
,
psi2
);
const
Real
omega_p
=
knorm
*
cp
;
const
Real
omega_s
=
knorm
*
cs
;
Real
phase_p
=
X
.
dot
(
kpush
)
-
omega_p
*
current_time
;
Real
phase_s
=
X
.
dot
(
kpush
)
-
omega_s
*
current_time
;
vel
=
A
*
(
kpush
*
omega_p
*
cos
(
phase_p
)
+
kshear1
*
omega_s
*
cos
(
phase_s
)
+
kshear2
*
omega_s
*
cos
(
phase_s
));
};
template
<
ElementType
_type
>
class
SolutionFunctor
:
public
BC
::
Dirichlet
::
DirichletFunctor
{
public
:
SolutionFunctor
(
Real
current_time
,
SolidMechanicsModel
&
model
)
:
current_time
(
current_time
),
model
(
model
)
{}
public
:
// static constexpr UInt dim = DimensionHelper<_type>::dim;
static
constexpr
UInt
dim
=
ElementClass
<
_type
>::
getSpatialDimension
();
inline
void
operator
()(
UInt
node
,
Vector
<
bool
>
&
flags
,
Vector
<
Real
>
&
primal
,
const
Vector
<
Real
>
&
coord
)
const
{
flags
(
0
)
=
true
;
auto
&
vel
=
model
.
getVelocity
();
auto
itVel
=
vel
.
begin
(
model
.
getSpatialDimension
());
Vector
<
Real
>
v
=
itVel
[
node
];
solution_disp
<
dim
>
(
primal
,
coord
,
current_time
);
solution_vel
<
dim
>
(
v
,
coord
,
current_time
);
}
private
:
Real
current_time
;
SolidMechanicsModel
&
model
;
};
template
<
ElementType
_type
,
typename
AM
>
void
test_body
(
SolidMechanicsModel
&
model
,
AM
analysis_method
)
{
// constexpr UInt dim = DimensionHelper<_type>::dim;
static
constexpr
UInt
dim
=
ElementClass
<
_type
>::
getSpatialDimension
();
getStaticParser
().
parse
(
"test_solid_mechanics_model_"
"dynamics_material.dat"
);
model
.
initFull
(
SolidMechanicsModelOptions
(
analysis_method
));
bool
dump_paraview
=
false
;
if
(
dump_paraview
)
{
std
::
stringstream
base_name
;
base_name
<<
"bar"
<<
analysis_method
<<
_type
;
model
.
setBaseName
(
base_name
.
str
());
model
.
addDumpFieldVector
(
"displacement"
);
model
.
addDumpField
(
"mass"
);
model
.
addDumpField
(
"velocity"
);
model
.
addDumpField
(
"acceleration"
);
model
.
addDumpFieldVector
(
"external_force"
);
model
.
addDumpFieldVector
(
"internal_force"
);
model
.
addDumpField
(
"stress"
);
model
.
addDumpField
(
"strain"
);
}
Real
time_step
=
model
.
getStableTimeStep
()
/
10.
;
model
.
setTimeStep
(
time_step
);
std
::
cout
<<
"timestep: "
<<
time_step
<<
std
::
endl
;
UInt
nb_nodes
=
model
.
getMesh
().
getNbNodes
();
UInt
spatial_dimension
=
model
.
getSpatialDimension
();
auto
&
nodes
=
model
.
getMesh
().
getNodes
();
auto
&
disp
=
model
.
getDisplacement
();
auto
&
vel
=
model
.
getVelocity
();
Array
<
Real
>
disp_solution
(
nb_nodes
,
spatial_dimension
);
Real
current_time
=
0
;
auto
itNodes
=
nodes
.
begin
(
spatial_dimension
);
auto
itDisp
=
disp
.
begin
(
spatial_dimension
);
auto
itVel
=
vel
.
begin
(
spatial_dimension
);
for
(
UInt
n
=
0
;
n
<
nb_nodes
;
++
n
,
++
itNodes
,
++
itDisp
,
++
itVel
)
{
solution_disp
<
dim
>
(
*
itDisp
,
*
itNodes
,
current_time
);
solution_vel
<
dim
>
(
*
itVel
,
*
itNodes
,
current_time
);
}
if
(
dump_paraview
)
model
.
dump
();
/// boundary conditions
model
.
applyBC
(
SolutionFunctor
<
_type
>
(
current_time
,
model
),
"BC"
);
Real
max_error
=
0.
;
Real
wave_velocity
=
1.
;
// sqrt(E/rho) = sqrt(1/1) = 1
Real
simulation_time
=
5
/
wave_velocity
;
UInt
max_steps
=
simulation_time
/
time_step
;
// 100
UInt
ndump
=
50
;
UInt
dump_freq
=
max_steps
/
ndump
;
std
::
cout
<<
"max_steps: "
<<
max_steps
<<
std
::
endl
;
for
(
UInt
s
=
0
;
s
<
max_steps
;
++
s
,
current_time
+=
time_step
)
{
if
(
s
%
dump_freq
==
0
&&
dump_paraview
)
model
.
dump
();
/// boundary conditions
model
.
applyBC
(
SolutionFunctor
<
_type
>
(
current_time
,
model
),
"BC"
);
// compute the disp solution
auto
itDispSolution
=
disp_solution
.
begin
(
spatial_dimension
);
itNodes
=
nodes
.
begin
(
spatial_dimension
);
for
(
UInt
n
=
0
;
n
<
nb_nodes
;
++
n
,
++
itNodes
,
++
itDispSolution
)
{
solution_disp
<
dim
>
(
*
itDispSolution
,
*
itNodes
,
current_time
);
}
// compute the error solution
itDispSolution
=
disp_solution
.
begin
(
spatial_dimension
);
itDisp
=
disp
.
begin
(
spatial_dimension
);
Real
disp_error
=
0.
;
for
(
UInt
n
=
0
;
n
<
nb_nodes
;
++
n
,
++
itDispSolution
,
++
itDisp
)
{
auto
diff
=
*
itDispSolution
-
*
itDisp
;
for
(
UInt
i
=
0
;
i
<
spatial_dimension
;
++
i
)
{
disp_error
+=
diff
(
i
)
*
diff
(
i
);
}
}
disp_error
=
sqrt
(
disp_error
)
/
nb_nodes
;
max_error
=
std
::
max
(
disp_error
,
max_error
);
ASSERT_NEAR
(
disp_error
,
0.
,
1e-3
);
model
.
solveStep
();
}
std
::
cout
<<
"max error: "
<<
max_error
<<
std
::
endl
;
}
#ifdef AKANTU_IMPLICIT
TYPED_TEST
(
TestSMMFixtureBar
,
DynamicsImplicit
)
{
test_body
<
this
->
type
>
(
*
(
this
->
model
),
_implicit_dynamic
);
}
#endif
TYPED_TEST
(
TestSMMFixtureBar
,
DynamicsExplicit
)
{
test_body
<
TestFixture
::
type
>
(
*
(
this
->
model
),
_explicit_lumped_mass
);
}
}
Event Timeline
Log In to Comment