Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F77910994
main.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
Sat, Aug 17, 03:56
Size
8 KB
Mime Type
text/x-c
Expires
Mon, Aug 19, 03:56 (2 d)
Engine
blob
Format
Raw Data
Handle
19946618
Attached To
rGOOSEFEM GooseFEM
main.cpp
View Options
#include <GMatElastoPlastic/Cartesian3d.h>
#include <GooseFEM/GooseFEM.h>
#include <GooseFEM/MatrixPartitionedTyings.h>
#include <GooseFEM/VectorPartitionedTyings.h>
#include <GooseFEM/TyingsPeriodic.h>
#include <XDMFWrite_HighFive.hpp>
#include <highfive/H5Easy.hpp>
namespace
GM
=
GMatElastoPlastic
::
Cartesian3d
;
namespace
GF
=
GooseFEM
;
namespace
PV
=
XDMFWrite_HighFive
;
namespace
H5
=
H5Easy
;
int
main
()
{
// mesh
// ----
// define mesh
GF
::
Mesh
::
Quad4
::
Regular
mesh
(
5
,
5
);
// mesh dimensions
size_t
nelem
=
mesh
.
nelem
();
size_t
nne
=
mesh
.
nne
();
size_t
ndim
=
mesh
.
ndim
();
// mesh definitions
xt
::
xtensor
<
double
,
2
>
coor
=
mesh
.
coor
();
xt
::
xtensor
<
size_t
,
2
>
conn
=
mesh
.
conn
();
xt
::
xtensor
<
size_t
,
2
>
dofs
=
mesh
.
dofs
();
xt
::
xtensor
<
size_t
,
2
>
elmat
=
mesh
.
elementgrid
();
// periodicity and fixed displacements DOFs
// ----------------------------------------
// add control nodes
GF
::
Tyings
::
Control
control
(
coor
,
dofs
);
coor
=
control
.
coor
();
dofs
=
control
.
dofs
();
xt
::
xtensor
<
size_t
,
2
>
control_dofs
=
control
.
controlDofs
();
xt
::
xtensor
<
size_t
,
1
>
control_nodes
=
control
.
controlNodes
();
// extract fixed DOFs:
// - all control nodes: to prescribe the deformation gradient
// - one node of the mesh: to remove rigid body modes
xt
::
xtensor
<
size_t
,
1
>
iip
=
xt
::
concatenate
(
xt
::
xtuple
(
xt
::
reshape_view
(
control_dofs
,
{
ndim
*
ndim
}),
xt
::
reshape_view
(
xt
::
view
(
dofs
,
xt
::
keep
(
mesh
.
nodesOrigin
()),
xt
::
all
()),
{
ndim
})));
// get DOF-tyings, reorganise system
GF
::
Tyings
::
Periodic
tyings
(
coor
,
dofs
,
control_dofs
,
mesh
.
nodesPeriodic
(),
iip
);
dofs
=
tyings
.
dofs
();
// simulation variables
// --------------------
// vector definition:
// provides methods to switch between dofval/nodeval/elemvec, or to manipulate a part of them
GF
::
VectorPartitionedTyings
vector
(
conn
,
dofs
,
tyings
.
Cdu
(),
tyings
.
Cdp
(),
tyings
.
Cdi
());
// nodal quantities
xt
::
xtensor
<
double
,
2
>
disp
=
xt
::
zeros
<
double
>
(
coor
.
shape
());
// nodal displacement
xt
::
xtensor
<
double
,
2
>
du
=
xt
::
zeros
<
double
>
(
coor
.
shape
());
// iterative displacement update
xt
::
xtensor
<
double
,
2
>
fint
=
xt
::
zeros
<
double
>
(
coor
.
shape
());
// internal force
xt
::
xtensor
<
double
,
2
>
fext
=
xt
::
zeros
<
double
>
(
coor
.
shape
());
// external force
xt
::
xtensor
<
double
,
2
>
fres
=
xt
::
zeros
<
double
>
(
coor
.
shape
());
// residual force
// element vectors / matrix
xt
::
xtensor
<
double
,
3
>
ue
=
xt
::
empty
<
double
>
({
nelem
,
nne
,
ndim
});
xt
::
xtensor
<
double
,
3
>
fe
=
xt
::
empty
<
double
>
({
nelem
,
nne
,
ndim
});
xt
::
xtensor
<
double
,
3
>
Ke
=
xt
::
empty
<
double
>
({
nelem
,
nne
*
ndim
,
nne
*
ndim
});
// DOF values
xt
::
xtensor
<
double
,
1
>
Fext
=
xt
::
zeros
<
double
>
({
tyings
.
nni
()});
xt
::
xtensor
<
double
,
1
>
Fint
=
xt
::
zeros
<
double
>
({
tyings
.
nni
()});
// element/material definition
// ---------------------------
// FEM quadrature
GF
::
Element
::
Quad4
::
QuadraturePlanar
elem
(
vector
.
AsElement
(
coor
));
size_t
nip
=
elem
.
nip
();
xt
::
xtensor
<
double
,
4
>
dV
=
elem
.
AsTensor
<
2
>
(
elem
.
dV
());
// material model
// even though the problem is 2-d, the material model is 3-d, plane strain is implicitly assumed
GM
::
Array
<
2
>
mat
({
nelem
,
nip
});
size_t
tdim
=
3
;
// some artificial material definition
xt
::
xtensor
<
size_t
,
1
>
ehard
=
xt
::
ravel
(
xt
::
view
(
elmat
,
xt
::
range
(
0
,
2
),
xt
::
range
(
0
,
2
)));
xt
::
xtensor
<
size_t
,
2
>
Ihard
=
xt
::
zeros
<
size_t
>
({
nelem
,
nip
});
xt
::
view
(
Ihard
,
xt
::
keep
(
ehard
),
xt
::
all
())
=
1ul
;
xt
::
xtensor
<
size_t
,
2
>
Isoft
=
xt
::
ones
<
size_t
>
({
nelem
,
nip
})
-
Ihard
;
mat
.
setLinearHardening
(
Isoft
,
1.0
,
1.0
,
0.05
,
0.05
);
mat
.
setElastic
(
Ihard
,
1.0
,
1.0
);
// solve
// -----
// allocate tensors
xt
::
xtensor
<
double
,
4
>
Eps
=
xt
::
empty
<
double
>
({
nelem
,
nip
,
tdim
,
tdim
});
xt
::
xtensor
<
double
,
4
>
Sig
=
xt
::
empty
<
double
>
({
nelem
,
nip
,
tdim
,
tdim
});
xt
::
xtensor
<
double
,
6
>
C
=
xt
::
empty
<
double
>
({
nelem
,
nip
,
tdim
,
tdim
,
tdim
,
tdim
});
// allocate system matrix
GF
::
MatrixPartitionedTyings
K
(
conn
,
dofs
,
tyings
.
Cdu
(),
tyings
.
Cdp
());
GF
::
MatrixPartitionedTyingsSolver
<>
Solver
;
// allocate internal variables
double
res
;
// some shear strain history
xt
::
xtensor
<
double
,
1
>
dgamma
=
0.001
*
xt
::
ones
<
double
>
({
101
});
dgamma
(
0
)
=
0.0
;
// initialise output
// - output-file containing data
HighFive
::
File
file
(
"main.h5"
,
HighFive
::
File
::
Overwrite
);
// - ParaView meta-data
PV
::
TimeSeries
xdmf
;
// - write mesh
H5
::
dump
(
file
,
"/conn"
,
conn
);
H5
::
dump
(
file
,
"/coor"
,
coor
);
// loop over increments
for
(
size_t
inc
=
0
;
inc
<
dgamma
.
size
();
++
inc
)
{
// update history
mat
.
increment
();
// iterate
for
(
size_t
iter
=
0
;;
++
iter
)
{
// strain
vector
.
asElement
(
disp
,
ue
);
elem
.
symGradN_vector
(
ue
,
Eps
);
// stress & tangent
mat
.
setStrain
(
Eps
);
mat
.
stress
(
Sig
);
mat
.
tangent
(
C
);
// internal force
elem
.
int_gradN_dot_tensor2_dV
(
Sig
,
fe
);
vector
.
assembleNode
(
fe
,
fint
);
// stiffness matrix
elem
.
int_gradN_dot_tensor4_dot_gradNT_dV
(
C
,
Ke
);
K
.
assemble
(
Ke
);
// residual
xt
::
noalias
(
fres
)
=
fext
-
fint
;
// check for convergence (skip the zeroth iteration, as the residual still vanishes)
if
(
iter
>
0
)
{
// - internal/external force as DOFs (account for periodicity)
vector
.
asDofs_i
(
fext
,
Fext
);
vector
.
asDofs_i
(
fint
,
Fint
);
// - extract reaction force
vector
.
copy_p
(
Fint
,
Fext
);
// - norm of the residual and the reaction force
double
nfres
=
xt
::
sum
(
xt
::
abs
(
Fext
-
Fint
))[
0
];
double
nfext
=
xt
::
sum
(
xt
::
abs
(
Fext
))[
0
];
// - relative residual, for convergence check
if
(
nfext
)
res
=
nfres
/
nfext
;
else
res
=
nfres
;
// - print progress to screen
std
::
cout
<<
"inc = "
<<
inc
<<
", "
<<
"iter = "
<<
iter
<<
", "
<<
"res = "
<<
res
<<
std
::
endl
;
// - check for convergence
if
(
res
<
1.0e-5
)
{
break
;
}
// - safe-guard from infinite loop
if
(
iter
>
20
)
{
throw
std
::
runtime_error
(
"Maximal number of iterations exceeded"
);
}
}
// initialise displacement update
du
.
fill
(
0.0
);
// set fixed displacements
if
(
iter
==
0
)
{
du
(
control_nodes
(
0
),
1
)
=
dgamma
(
inc
);
}
// solve
Solver
.
solve
(
K
,
fres
,
du
);
// add displacement update
disp
+=
du
;
}
// post-process
// - compute strain and stress
vector
.
asElement
(
disp
,
ue
);
elem
.
symGradN_vector
(
ue
,
Eps
);
mat
.
setStrain
(
Eps
);
mat
.
stress
(
Sig
);
// - element average stress
xt
::
xtensor
<
double
,
3
>
Sigelem
=
xt
::
average
(
Sig
,
dV
,
{
1
});
xt
::
xtensor
<
double
,
3
>
Epselem
=
xt
::
average
(
Eps
,
dV
,
{
1
});
// - macroscopic strain and stress
xt
::
xtensor_fixed
<
double
,
xt
::
xshape
<
3
,
3
>>
Sigbar
=
xt
::
average
(
Sig
,
dV
,
{
0
,
1
});
xt
::
xtensor_fixed
<
double
,
xt
::
xshape
<
3
,
3
>>
Epsbar
=
xt
::
average
(
Eps
,
dV
,
{
0
,
1
});
// - write to output-file: increment numbers
H5
::
dump
(
file
,
"/stored"
,
inc
,
{
inc
});
// - write to output-file: macroscopic response
H5
::
dump
(
file
,
"/macroscopic/sigeq"
,
GM
::
Sigeq
(
Sigbar
)(),
{
inc
});
H5
::
dump
(
file
,
"/macroscopic/epseq"
,
GM
::
Epseq
(
Epsbar
)(),
{
inc
});
// - write to output-file: element quantities
H5
::
dump
(
file
,
"/sigeq/"
+
std
::
to_string
(
inc
),
GM
::
Sigeq
(
Sigelem
)());
H5
::
dump
(
file
,
"/epseq/"
+
std
::
to_string
(
inc
),
GM
::
Epseq
(
Epselem
)());
H5
::
dump
(
file
,
"/disp/"
+
std
::
to_string
(
inc
),
GF
::
as3d
(
disp
));
// - update ParaView meta-data
// xdmf.push_back({
// PV::Topology(file, "/conn", mesh.getElementType()),
// PV::Geometry(file, "/coor"),
// PV::Attribute(file, "/disp/" + std::to_string(inc), PV::AttributeCenter::Node, "Displacement"),
// PV::Attribute(file, "/sigeq/" + std::to_string(inc), PV::AttributeCenter::Cell, "Eq. stress"),
// PV::Attribute(file, "/epseq/" + std::to_string(inc), PV::AttributeCenter::Cell, "Eq. strain")});
}
// write ParaView meta-data
PV
::
write
(
"main.xdmf"
,
xdmf
.
get
());
return
0
;
}
Event Timeline
Log In to Comment