Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F73596404
patch_test_explicit_anisotropic.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
Tue, Jul 23, 00:31
Size
10 KB
Mime Type
text/x-c
Expires
Thu, Jul 25, 00:31 (2 d)
Engine
blob
Format
Raw Data
Handle
19231884
Attached To
rAKA akantu
patch_test_explicit_anisotropic.cc
View Options
/**
* @file patch_test_explicit.cc
*
* @author David Simon Kammer <david.kammer@epfl.ch>
* @author Nicolas Richart <nicolas.richart@epfl.ch>
* @author Cyprien Wolff <cyprien.wolff@epfl.ch>
*
* @date Thu Feb 17 16:05:48 2011
*
* @brief patch test for elastic material in solid mechanics model
*
* @section LICENSE
*
* Copyright (©) 2010-2011 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 <iostream>
#include "solid_mechanics_model.hh"
using
namespace
akantu
;
Real
alpha
[
3
][
4
]
=
{
{
0.01
,
0.02
,
0.03
,
0.04
},
{
0.05
,
0.06
,
0.07
,
0.08
},
{
0.09
,
0.10
,
0.11
,
0.12
}
};
//Stiffness tensor, rotated by hand
Real
C
[
3
][
3
][
3
][
3
]
=
{{{{
112.93753505
,
1.85842452538e-10
,
-
4.47654358027e-10
},
{
1.85847317471e-10
,
54.2334345331
,
-
3.69840984824
},
{
-
4.4764768395e-10
,
-
3.69840984824
,
56.848605217
}},
{{
1.85847781609e-10
,
25.429294233
,
-
3.69840984816
},
{
25.429294233
,
3.31613847493e-10
,
-
8.38797920011e-11
},
{
-
3.69840984816
,
-
8.38804581349e-11
,
-
1.97875715813e-10
}},
{{
-
4.47654358027e-10
,
-
3.69840984816
,
28.044464917
},
{
-
3.69840984816
,
2.09374961813e-10
,
9.4857455224e-12
},
{
28.044464917
,
9.48308098714e-12
,
-
2.1367885239e-10
}}},
{{{
1.85847781609e-10
,
25.429294233
,
-
3.69840984816
},
{
25.429294233
,
3.31613847493e-10
,
-
8.38793479119e-11
},
{
-
3.69840984816
,
-
8.38795699565e-11
,
-
1.97876381947e-10
}},
{{
54.2334345331
,
3.31617400207e-10
,
2.09372075233e-10
},
{
3.3161562385e-10
,
115.552705733
,
-
3.15093728886e-10
},
{
2.09372075233e-10
,
-
3.15090176173e-10
,
54.2334345333
}},
{{
-
3.69840984824
,
-
8.38795699565e-11
,
9.48219280872e-12
},
{
-
8.38795699565e-11
,
-
3.1509195253e-10
,
25.4292942335
},
{
9.48441325477e-12
,
25.4292942335
,
3.69840984851
}}},
{{{
-
4.47653469848e-10
,
-
3.69840984816
,
28.044464917
},
{
-
3.69840984816
,
2.09374073634e-10
,
9.48752187924e-12
},
{
28.044464917
,
9.48552347779e-12
,
-
2.1367885239e-10
}},
{{
-
3.69840984824
,
-
8.3884899027e-11
,
9.48219280872e-12
},
{
-
8.3884899027e-11
,
-
3.150972816e-10
,
25.4292942335
},
{
9.48041645188e-12
,
25.4292942335
,
3.69840984851
}},
{{
56.848605217
,
-
1.97875493768e-10
,
-
2.13681516925e-10
},
{
-
1.97877270125e-10
,
54.2334345333
,
3.69840984851
},
{
-
2.13683293282e-10
,
3.69840984851
,
112.93753505
}}}};
/* -------------------------------------------------------------------------- */
template
<
ElementType
type
>
static
Matrix
<
Real
>
prescribed_grad_u
()
{
UInt
spatial_dimension
=
ElementClass
<
type
>::
getSpatialDimension
();
Matrix
<
Real
>
grad_u
(
spatial_dimension
,
spatial_dimension
);
for
(
UInt
i
=
0
;
i
<
spatial_dimension
;
++
i
)
{
for
(
UInt
j
=
0
;
j
<
spatial_dimension
;
++
j
)
{
grad_u
(
i
,
j
)
=
alpha
[
i
][
j
+
1
];
}
}
return
grad_u
;
}
template
<
ElementType
type
>
static
Matrix
<
Real
>
prescribed_stress
()
{
UInt
spatial_dimension
=
ElementClass
<
type
>::
getSpatialDimension
();
Matrix
<
Real
>
stress
(
spatial_dimension
,
spatial_dimension
);
//plane strain in 2d
Matrix
<
Real
>
strain
(
spatial_dimension
,
spatial_dimension
);
Matrix
<
Real
>
pstrain
;
pstrain
=
prescribed_grad_u
<
type
>
();
Real
trace
=
0
;
/// symetric part of the strain tensor
for
(
UInt
i
=
0
;
i
<
spatial_dimension
;
++
i
)
for
(
UInt
j
=
0
;
j
<
spatial_dimension
;
++
j
)
strain
(
i
,
j
)
=
0.5
*
(
pstrain
(
i
,
j
)
+
pstrain
(
j
,
i
));
for
(
UInt
i
=
0
;
i
<
spatial_dimension
;
++
i
)
trace
+=
strain
(
i
,
i
);
for
(
UInt
i
=
0
;
i
<
spatial_dimension
;
++
i
)
{
for
(
UInt
j
=
0
;
j
<
spatial_dimension
;
++
j
)
{
stress
(
i
,
j
)
=
0
;
for
(
UInt
k
=
0
;
k
<
spatial_dimension
;
++
k
)
{
for
(
UInt
l
=
0
;
l
<
spatial_dimension
;
++
l
)
{
stress
(
i
,
j
)
+=
C
[
i
][
j
][
k
][
l
]
*
strain
(
k
,
l
);
}
}
}
}
return
stress
;
}
#define TYPE _tetrahedron_4
/* -------------------------------------------------------------------------- */
int
main
(
int
argc
,
char
*
argv
[])
{
initialize
(
"material_anisotropic.dat"
,
argc
,
argv
);
UInt
dim
=
ElementClass
<
TYPE
>::
getSpatialDimension
();
const
ElementType
element_type
=
TYPE
;
UInt
damping_steps
=
600000
;
UInt
damping_interval
=
50
;
Real
damping_ratio
=
0.99
;
UInt
additional_steps
=
20000
;
UInt
max_steps
=
damping_steps
+
additional_steps
;
/// load mesh
Mesh
my_mesh
(
dim
);
std
::
stringstream
filename
;
filename
<<
TYPE
<<
".msh"
;
my_mesh
.
read
(
filename
.
str
());
UInt
nb_nodes
=
my_mesh
.
getNbNodes
();
/// declaration of model
SolidMechanicsModel
my_model
(
my_mesh
);
/// model initialization
my_model
.
initFull
(
SolidMechanicsModelOptions
(
_explicit_lumped_mass
));
std
::
cout
<<
my_model
.
getMaterial
(
0
)
<<
std
::
endl
;
Real
time_step
=
my_model
.
getStableTimeStep
()
/
5.
;
my_model
.
setTimeStep
(
time_step
);
my_model
.
assembleMassLumped
();
std
::
cout
<<
"The number of time steps is: "
<<
max_steps
<<
" ("
<<
time_step
<<
"s)"
<<
std
::
endl
;
// boundary conditions
const
Array
<
Real
>
&
coordinates
=
my_mesh
.
getNodes
();
Array
<
Real
>
&
displacement
=
my_model
.
getDisplacement
();
Array
<
bool
>
&
boundary
=
my_model
.
getBlockedDOFs
();
MeshUtils
::
buildFacets
(
my_mesh
);
my_mesh
.
createBoundaryGroupFromGeometry
();
// Loop over (Sub)Boundar(ies)
// Loop over (Sub)Boundar(ies)
for
(
GroupManager
::
const_element_group_iterator
it
(
my_mesh
.
element_group_begin
());
it
!=
my_mesh
.
element_group_end
();
++
it
)
{
for
(
ElementGroup
::
const_node_iterator
nodes_it
(
it
->
second
->
node_begin
());
nodes_it
!=
it
->
second
->
node_end
();
++
nodes_it
)
{
UInt
n
(
*
nodes_it
);
std
::
cout
<<
"Node "
<<
*
nodes_it
<<
std
::
endl
;
for
(
UInt
i
=
0
;
i
<
dim
;
++
i
)
{
displacement
(
n
,
i
)
=
alpha
[
i
][
0
];
for
(
UInt
j
=
0
;
j
<
dim
;
++
j
)
{
displacement
(
n
,
i
)
+=
alpha
[
i
][
j
+
1
]
*
coordinates
(
n
,
j
);
}
boundary
(
n
,
i
)
=
true
;
}
}
}
// Actually, loop over all nodes, since I wanna test a static solution
for
(
UInt
n
=
0
;
n
<
nb_nodes
;
++
n
)
{
for
(
UInt
i
=
0
;
i
<
dim
;
++
i
)
{
displacement
(
n
,
i
)
=
alpha
[
i
][
0
];
for
(
UInt
j
=
0
;
j
<
dim
;
++
j
)
{
displacement
(
n
,
i
)
+=
alpha
[
i
][
j
+
1
]
*
coordinates
(
n
,
j
);
}
}
}
Array
<
Real
>
&
velocity
=
my_model
.
getVelocity
();
std
::
ofstream
energy
;
std
::
stringstream
energy_filename
;
energy_filename
<<
"energy_"
<<
TYPE
<<
".csv"
;
energy
.
open
(
energy_filename
.
str
().
c_str
());
energy
<<
"id,time,ekin"
<<
std
::
endl
;
Real
ekin_mean
=
0.
;
/* ------------------------------------------------------------------------ */
/* Main loop */
/* ------------------------------------------------------------------------ */
UInt
s
;
for
(
s
=
1
;
s
<=
max_steps
;
++
s
)
{
if
(
s
%
10000
==
0
)
std
::
cout
<<
"passing step "
<<
s
<<
"/"
<<
max_steps
<<
" ("
<<
s
*
time_step
<<
"s)"
<<
std
::
endl
;
// damp velocity in order to find equilibrium
if
((
s
<
damping_steps
)
&&
(
s
%
damping_interval
==
0
))
{
velocity
*=
damping_ratio
;
}
if
(
s
%
1000
==
0
)
{
ekin_mean
=
ekin_mean
/
1000.
;
std
::
cout
<<
"Ekin mean = "
<<
ekin_mean
<<
std
::
endl
;
if
(
ekin_mean
<
1e-10
)
break
;
ekin_mean
=
0.
;
}
my_model
.
explicitPred
();
my_model
.
updateResidual
();
my_model
.
updateAcceleration
();
my_model
.
explicitCorr
();
akantu
::
Real
ekin
=
my_model
.
getKineticEnergy
();
ekin_mean
+=
ekin
;
if
(
s
%
1000
==
0
)
energy
<<
s
<<
","
<<
s
*
time_step
<<
","
<<
ekin
<<
std
::
endl
;
}
energy
.
close
();
UInt
nb_quadrature_points
=
my_model
.
getFEEngine
().
getNbQuadraturePoints
(
TYPE
);
Array
<
Real
>
&
stress_vect
=
const_cast
<
Array
<
Real
>
&>
(
my_model
.
getMaterial
(
0
).
getStress
(
element_type
));
Array
<
Real
>
&
gradu_vect
=
const_cast
<
Array
<
Real
>
&>
(
my_model
.
getMaterial
(
0
).
getGradU
(
element_type
));
Array
<
Real
>::
iterator
<
Matrix
<
Real
>
>
stress_it
=
stress_vect
.
begin
(
dim
,
dim
);
Array
<
Real
>::
iterator
<
Matrix
<
Real
>
>
gradu_it
=
gradu_vect
.
begin
(
dim
,
dim
);
Matrix
<
Real
>
presc_stress
;
presc_stress
=
prescribed_stress
<
TYPE
>
();
Matrix
<
Real
>
presc_gradu
;
presc_gradu
=
prescribed_grad_u
<
TYPE
>
();
UInt
nb_element
=
my_mesh
.
getNbElement
(
TYPE
);
Real
gradu_tolerance
=
1e-9
;
Real
stress_tolerance
=
1e-2
;
if
(
s
>
max_steps
)
{
stress_tolerance
=
1e-4
;
gradu_tolerance
=
1e-7
;
}
for
(
UInt
el
=
0
;
el
<
nb_element
;
++
el
)
{
for
(
UInt
q
=
0
;
q
<
nb_quadrature_points
;
++
q
)
{
Matrix
<
Real
>
&
stress
=
*
stress_it
;
Matrix
<
Real
>
&
gradu
=
*
gradu_it
;
for
(
UInt
i
=
0
;
i
<
dim
;
++
i
)
{
for
(
UInt
j
=
0
;
j
<
dim
;
++
j
)
{
if
(
!
(
std
::
abs
(
gradu
(
i
,
j
)
-
presc_gradu
(
i
,
j
))
<
gradu_tolerance
))
{
std
::
cerr
<<
"gradu["
<<
i
<<
","
<<
j
<<
"] = "
<<
gradu
(
i
,
j
)
<<
" but should be = "
<<
presc_gradu
(
i
,
j
)
<<
" (-"
<<
std
::
abs
(
gradu
(
i
,
j
)
-
presc_gradu
(
i
,
j
))
<<
") [el : "
<<
el
<<
" - q : "
<<
q
<<
"]"
<<
std
::
endl
;
std
::
cerr
<<
gradu
<<
presc_gradu
<<
std
::
endl
;
return
EXIT_FAILURE
;
}
if
(
!
(
std
::
abs
(
stress
(
i
,
j
)
-
presc_stress
(
i
,
j
))
<
stress_tolerance
))
{
std
::
cerr
<<
"stress["
<<
i
<<
","
<<
j
<<
"] = "
<<
stress
(
i
,
j
)
<<
" but should be = "
<<
presc_stress
(
i
,
j
)
<<
" (-"
<<
std
::
abs
(
stress
(
i
,
j
)
-
presc_stress
(
i
,
j
))
<<
") [el : "
<<
el
<<
" - q : "
<<
q
<<
"]"
<<
std
::
endl
;
std
::
cerr
<<
stress
<<
presc_stress
<<
std
::
endl
;
return
EXIT_FAILURE
;
}
}
}
++
stress_it
;
++
gradu_it
;
}
}
for
(
UInt
n
=
0
;
n
<
nb_nodes
;
++
n
)
{
for
(
UInt
i
=
0
;
i
<
dim
;
++
i
)
{
Real
disp
=
alpha
[
i
][
0
];
for
(
UInt
j
=
0
;
j
<
dim
;
++
j
)
{
disp
+=
alpha
[
i
][
j
+
1
]
*
coordinates
(
n
,
j
);
}
if
(
!
(
std
::
abs
(
displacement
(
n
,
i
)
-
disp
)
<
1e-7
))
{
std
::
cerr
<<
"displacement("
<<
n
<<
", "
<<
i
<<
")="
<<
displacement
(
n
,
i
)
<<
" should be equal to "
<<
disp
<<
"("
<<
displacement
(
n
,
i
)
-
disp
<<
")"
<<
std
::
endl
;
return
EXIT_FAILURE
;
}
}
}
finalize
();
return
EXIT_SUCCESS
;
}
Event Timeline
Log In to Comment