Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F74569757
main_velocityVerlet.cpp
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
Sun, Jul 28, 12:20
Size
4 KB
Mime Type
text/x-c
Expires
Tue, Jul 30, 12:20 (2 d)
Engine
blob
Format
Raw Data
Handle
19411805
Attached To
rGOOSEFEM GooseFEM
main_velocityVerlet.cpp
View Options
#include <GooseFEM/GooseFEM.h>
#include <GMatElastoPlasticQPot/Cartesian2d.h>
#include <LowFive.h>
// -------------------------------------------------------------------------------------------------
namespace
GF
=
GooseFEM
;
namespace
QD
=
GooseFEM
::
Element
::
Quad4
;
namespace
GM
=
GMatElastoPlasticQPot
::
Cartesian2d
;
// -------------------------------------------------------------------------------------------------
inline
double
sqdot
(
const
xt
::
xtensor
<
double
,
1
>
&
M
,
const
xt
::
xtensor
<
double
,
1
>
&
V
)
{
double
out
=
0.
;
for
(
size_t
i
=
0
;
i
<
M
.
size
()
;
++
i
)
out
+=
M
(
i
)
*
V
(
i
)
*
V
(
i
);
return
out
;
}
// -------------------------------------------------------------------------------------------------
int
main
()
{
// simulation parameters
double
T
=
60.
;
// total time
double
dt
=
1.e-2
;
// time increment
size_t
nx
=
60
;
// number of elements in both directions
double
gamma
=
.05
;
// displacement step
// get mesh & quadrature
GF
::
Mesh
::
Quad4
::
Regular
mesh
(
nx
,
nx
,
1.
);
xt
::
xtensor
<
double
,
2
>
coor
=
mesh
.
coor
();
xt
::
xtensor
<
size_t
,
2
>
conn
=
mesh
.
conn
();
xt
::
xtensor
<
size_t
,
2
>
dofs
=
mesh
.
dofsPeriodic
();
GF
::
Vector
vector
(
conn
,
dofs
);
QD
::
Quadrature
quad
(
vector
.
asElement
(
coor
));
size_t
nnode
=
mesh
.
nnode
();
size_t
ndim
=
mesh
.
ndim
();
size_t
nelem
=
mesh
.
nelem
();
size_t
nip
=
quad
.
nip
();
xt
::
xtensor
<
double
,
2
>
u
=
xt
::
zeros
<
double
>
(
coor
.
shape
());
xt
::
xtensor
<
double
,
2
>
v
=
xt
::
zeros
<
double
>
(
coor
.
shape
());
xt
::
xtensor
<
double
,
2
>
a
=
xt
::
zeros
<
double
>
(
coor
.
shape
());
xt
::
xtensor
<
double
,
2
>
v_n
=
xt
::
zeros
<
double
>
(
coor
.
shape
());
xt
::
xtensor
<
double
,
2
>
a_n
=
xt
::
zeros
<
double
>
(
coor
.
shape
());
xt
::
xtensor
<
double
,
2
>
fint
=
xt
::
zeros
<
double
>
(
coor
.
shape
());
xt
::
xtensor
<
double
,
2
>
fext
=
xt
::
zeros
<
double
>
(
coor
.
shape
());
xt
::
xtensor
<
double
,
2
>
fres
=
xt
::
zeros
<
double
>
(
coor
.
shape
());
xt
::
xtensor
<
double
,
4
>
Eps
=
xt
::
zeros
<
double
>
({
nelem
,
nip
,
ndim
,
ndim
});
xt
::
xtensor
<
double
,
4
>
Sig
=
xt
::
zeros
<
double
>
({
nelem
,
nip
,
ndim
,
ndim
});
// material definition
GM
::
Matrix
material
({
nelem
,
nip
});
xt
::
xtensor
<
size_t
,
2
>
Ihard
=
xt
::
zeros
<
size_t
>
({
nelem
,
nip
});
xt
::
xtensor
<
size_t
,
2
>
Isoft
=
xt
::
ones
<
size_t
>
({
nelem
,
nip
});
for
(
size_t
e
=
0
;
e
<
nx
*
nx
/
4
;
++
e
)
for
(
size_t
q
=
0
;
q
<
nip
;
++
q
)
Ihard
(
e
,
q
)
=
1
;
Isoft
-=
Ihard
;
material
.
setElastic
(
Ihard
,
100.
,
10.
);
material
.
setElastic
(
Isoft
,
100.
,
1.
);
material
.
check
();
// mass matrix
QD
::
Quadrature
nodalQuad
(
vector
.
asElement
(
coor
),
QD
::
Nodal
::
xi
(),
QD
::
Nodal
::
w
());
xt
::
xtensor
<
double
,
2
>
rho
=
1.0
*
xt
::
ones
<
double
>
({
nelem
,
nip
});
GF
::
MatrixDiagonal
M
(
conn
,
dofs
);
M
.
assemble
(
nodalQuad
.
int_N_scalar_NT_dV
(
rho
)
);
xt
::
xtensor
<
double
,
1
>
mass
=
M
.
asDiagonal
();
// update in macroscopic deformation gradient
xt
::
xtensor
<
double
,
2
>
dFbar
=
xt
::
zeros
<
double
>
({
2
,
2
});
dFbar
(
0
,
1
)
=
gamma
;
for
(
size_t
n
=
0
;
n
<
nnode
;
++
n
)
for
(
size_t
j
=
0
;
j
<
ndim
;
++
j
)
for
(
size_t
k
=
0
;
k
<
ndim
;
++
k
)
u
(
n
,
j
)
+=
dFbar
(
j
,
k
)
*
(
coor
(
n
,
k
)
-
coor
(
0
,
k
)
);
// output variables
xt
::
xtensor
<
double
,
1
>
Epot
=
xt
::
zeros
<
double
>
({
static_cast
<
size_t
>
(
T
/
dt
)});
xt
::
xtensor
<
double
,
1
>
Ekin
=
xt
::
zeros
<
double
>
({
static_cast
<
size_t
>
(
T
/
dt
)});
xt
::
xtensor
<
double
,
1
>
t
=
xt
::
zeros
<
double
>
({
static_cast
<
size_t
>
(
T
/
dt
)});
// loop over increments
for
(
size_t
inc
=
0
;
inc
<
static_cast
<
size_t
>
(
Epot
.
size
())
;
++
inc
)
{
// store history
xt
::
noalias
(
v_n
)
=
v
;
xt
::
noalias
(
a_n
)
=
a
;
// new displacement
xt
::
noalias
(
u
)
=
u
+
dt
*
v
+
0.5
*
std
::
pow
(
dt
,
2.
)
*
a
;
// compute strain/strain, and corresponding internal
quad
.
symGradN_vector
(
vector
.
asElement
(
u
),
Eps
);
material
.
Sig
(
Eps
,
Sig
);
vector
.
assembleNode
(
quad
.
int_gradN_dot_tensor2_dV
(
Sig
),
fint
);
// estimate new velocity
xt
::
noalias
(
v
)
=
v_n
+
dt
*
a_n
;
// compute residual force & solve
xt
::
noalias
(
fres
)
=
fext
-
fint
;
M
.
solve
(
fres
,
a
);
// re-estimate new velocity
xt
::
noalias
(
v
)
=
v_n
+
.5
*
dt
*
(
a_n
+
a
);
// compute residual force & solve
xt
::
noalias
(
fres
)
=
fext
-
fint
;
M
.
solve
(
fres
,
a
);
// new velocity
xt
::
noalias
(
v
)
=
v_n
+
.5
*
dt
*
(
a_n
+
a
);
// compute residual force & solve
xt
::
noalias
(
fres
)
=
fext
-
fint
;
M
.
solve
(
fres
,
a
);
// store output variables
xt
::
xtensor
<
double
,
2
>
E
=
material
.
energy
(
Eps
);
xt
::
xtensor
<
double
,
2
>
dV
=
quad
.
dV
();
xt
::
xtensor
<
double
,
1
>
V
=
vector
.
asDofs
(
v
);
t
(
inc
)
=
static_cast
<
double
>
(
inc
)
*
dt
;
Ekin
(
inc
)
=
0.5
*
sqdot
(
mass
,
V
);
Epot
(
inc
)
=
xt
::
sum
(
E
*
dV
)[
0
];
}
// write output variables to file
HighFive
::
File
file
(
"example.hdf5"
,
HighFive
::
File
::
Overwrite
);
LowFive
::
xtensor
::
dump
(
file
,
"/global/Epot"
,
Epot
);
LowFive
::
xtensor
::
dump
(
file
,
"/global/Ekin"
,
Ekin
);
LowFive
::
xtensor
::
dump
(
file
,
"/global/t"
,
t
);
LowFive
::
xtensor
::
dump
(
file
,
"/mesh/conn"
,
conn
);
LowFive
::
xtensor
::
dump
(
file
,
"/mesh/coor"
,
coor
);
LowFive
::
xtensor
::
dump
(
file
,
"/mesh/disp"
,
u
);
return
0
;
}
Event Timeline
Log In to Comment