Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F73204676
test_elastic_materials.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
Fri, Jul 19, 14:06
Size
28 KB
Mime Type
text/x-c++
Expires
Sun, Jul 21, 14:06 (2 d)
Engine
blob
Format
Raw Data
Handle
19162010
Attached To
rAKA akantu
test_elastic_materials.cc
View Options
/* -------------------------------------------------------------------------- */
#include "material_elastic.hh"
#include "material_elastic_orthotropic.hh"
#include "solid_mechanics_model.hh"
#include "test_material_fixtures.hh"
/* -------------------------------------------------------------------------- */
#include <gtest/gtest.h>
#include <type_traits>
/* -------------------------------------------------------------------------- */
using
namespace
akantu
;
using
types
=
::
testing
::
Types
<
Traits
<
MaterialElastic
,
1
>
,
Traits
<
MaterialElastic
,
2
>
,
Traits
<
MaterialElastic
,
3
>
,
Traits
<
MaterialElasticOrthotropic
,
1
>
,
Traits
<
MaterialElasticOrthotropic
,
2
>
,
Traits
<
MaterialElasticOrthotropic
,
3
>
,
Traits
<
MaterialElasticLinearAnisotropic
,
1
>
,
Traits
<
MaterialElasticLinearAnisotropic
,
2
>
,
Traits
<
MaterialElasticLinearAnisotropic
,
3
>>
;
/* -------------------------------------------------------------------------- */
template
<>
void
FriendMaterial
<
MaterialElastic
<
1
>>::
testComputeStress
()
{
Real
E
=
3.
;
setParam
(
"E"
,
E
);
Matrix
<
Real
>
eps
=
{{
2
}};
Matrix
<
Real
>
sigma
(
1
,
1
);
Real
sigma_th
=
2
;
this
->
computeStressOnQuad
(
eps
,
sigma
,
sigma_th
);
auto
solution
=
E
*
eps
(
0
,
0
)
+
sigma_th
;
ASSERT_NEAR
(
sigma
(
0
,
0
),
solution
,
1e-14
);
}
/* -------------------------------------------------------------------------- */
template
<>
void
FriendMaterial
<
MaterialElastic
<
1
>>::
testEnergyDensity
()
{
Real
eps
=
2
,
sigma
=
2
;
Real
epot
=
0
;
this
->
computePotentialEnergyOnQuad
({{
eps
}},
{{
sigma
}},
epot
);
Real
solution
=
2
;
ASSERT_NEAR
(
epot
,
solution
,
1e-14
);
}
/* -------------------------------------------------------------------------- */
template
<>
void
FriendMaterial
<
MaterialElastic
<
1
>>::
testComputeTangentModuli
()
{
Real
E
=
2
;
setParam
(
"E"
,
E
);
Matrix
<
Real
>
tangent
(
1
,
1
);
this
->
computeTangentModuliOnQuad
(
tangent
);
ASSERT_NEAR
(
tangent
(
0
,
0
),
E
,
1e-14
);
}
/* -------------------------------------------------------------------------- */
template
<>
void
FriendMaterial
<
MaterialElastic
<
1
>>::
testCelerity
()
{
Real
E
=
3.
,
rho
=
2.
;
setParam
(
"E"
,
E
);
setParam
(
"rho"
,
rho
);
auto
wave_speed
=
this
->
getCelerity
(
Element
());
auto
solution
=
std
::
sqrt
(
E
/
rho
);
ASSERT_NEAR
(
wave_speed
,
solution
,
1e-14
);
}
/* -------------------------------------------------------------------------- */
template
<>
void
FriendMaterial
<
MaterialElastic
<
2
>>::
testComputeStress
()
{
Real
E
=
1.
;
Real
nu
=
.3
;
Real
sigma_th
=
0.3
;
// thermal stress
setParam
(
"E"
,
E
);
setParam
(
"nu"
,
nu
);
Matrix
<
Real
>
rotation_matrix
=
getRandomRotation2d
();
auto
grad_u
=
this
->
getDeviatoricStrain
(
1.
).
block
(
0
,
0
,
2
,
2
);
auto
grad_u_rot
=
this
->
applyRotation
(
grad_u
,
rotation_matrix
);
Matrix
<
Real
>
sigma_rot
(
2
,
2
);
this
->
computeStressOnQuad
(
grad_u_rot
,
sigma_rot
,
sigma_th
);
auto
sigma
=
this
->
reverseRotation
(
sigma_rot
,
rotation_matrix
);
Matrix
<
Real
>
identity
(
2
,
2
);
identity
.
eye
();
Matrix
<
Real
>
sigma_expected
=
0.5
*
E
/
(
1
+
nu
)
*
(
grad_u
+
grad_u
.
transpose
())
+
sigma_th
*
identity
;
auto
diff
=
sigma
-
sigma_expected
;
Real
stress_error
=
diff
.
norm
<
L_inf
>
();
ASSERT_NEAR
(
stress_error
,
0.
,
1e-14
);
}
/* -------------------------------------------------------------------------- */
template
<>
void
FriendMaterial
<
MaterialElastic
<
2
>>::
testEnergyDensity
()
{
Matrix
<
Real
>
sigma
=
{{
1
,
2
},
{
2
,
4
}};
Matrix
<
Real
>
eps
=
{{
1
,
0
},
{
0
,
1
}};
Real
epot
=
0
;
Real
solution
=
2.5
;
this
->
computePotentialEnergyOnQuad
(
eps
,
sigma
,
epot
);
ASSERT_NEAR
(
epot
,
solution
,
1e-14
);
}
/* -------------------------------------------------------------------------- */
template
<>
void
FriendMaterial
<
MaterialElastic
<
2
>>::
testComputeTangentModuli
()
{
Real
E
=
1.
;
Real
nu
=
.3
;
setParam
(
"E"
,
E
);
setParam
(
"nu"
,
nu
);
Matrix
<
Real
>
tangent
(
3
,
3
);
/* Plane Strain */
// clang-format off
Matrix
<
Real
>
solution
=
{
{
1
-
nu
,
nu
,
0
},
{
nu
,
1
-
nu
,
0
},
{
0
,
0
,
(
1
-
2
*
nu
)
/
2
},
};
// clang-format on
solution
*=
E
/
((
1
+
nu
)
*
(
1
-
2
*
nu
));
this
->
computeTangentModuliOnQuad
(
tangent
);
Real
tangent_error
=
(
tangent
-
solution
).
norm
<
L_2
>
();
ASSERT_NEAR
(
tangent_error
,
0
,
1e-14
);
/* Plane Stress */
this
->
plane_stress
=
true
;
this
->
updateInternalParameters
();
// clang-format off
solution
=
{
{
1
,
nu
,
0
},
{
nu
,
1
,
0
},
{
0
,
0
,
(
1
-
nu
)
/
2
},
};
// clang-format on
solution
*=
E
/
(
1
-
nu
*
nu
);
this
->
computeTangentModuliOnQuad
(
tangent
);
tangent_error
=
(
tangent
-
solution
).
norm
<
L_2
>
();
ASSERT_NEAR
(
tangent_error
,
0
,
1e-14
);
}
/* -------------------------------------------------------------------------- */
template
<>
void
FriendMaterial
<
MaterialElastic
<
2
>>::
testCelerity
()
{
Real
E
=
1.
;
Real
nu
=
.3
;
Real
rho
=
2
;
setParam
(
"E"
,
E
);
setParam
(
"nu"
,
nu
);
setParam
(
"rho"
,
rho
);
auto
push_wave_speed
=
this
->
getPushWaveSpeed
(
Element
());
auto
celerity
=
this
->
getCelerity
(
Element
());
Real
K
=
E
/
(
3
*
(
1
-
2
*
nu
));
Real
mu
=
E
/
(
2
*
(
1
+
nu
));
Real
sol
=
std
::
sqrt
((
K
+
4.
/
3
*
mu
)
/
rho
);
ASSERT_NEAR
(
push_wave_speed
,
sol
,
1e-14
);
ASSERT_NEAR
(
celerity
,
sol
,
1e-14
);
auto
shear_wave_speed
=
this
->
getShearWaveSpeed
(
Element
());
sol
=
std
::
sqrt
(
mu
/
rho
);
ASSERT_NEAR
(
shear_wave_speed
,
sol
,
1e-14
);
}
/* -------------------------------------------------------------------------- */
template
<>
void
FriendMaterial
<
MaterialElastic
<
3
>>::
testComputeStress
()
{
Real
E
=
1.
;
Real
nu
=
.3
;
Real
sigma_th
=
0.3
;
// thermal stress
setParam
(
"E"
,
E
);
setParam
(
"nu"
,
nu
);
Matrix
<
Real
>
rotation_matrix
=
getRandomRotation3d
();
auto
grad_u
=
this
->
getDeviatoricStrain
(
1.
);
auto
grad_u_rot
=
this
->
applyRotation
(
grad_u
,
rotation_matrix
);
Matrix
<
Real
>
sigma_rot
(
3
,
3
);
this
->
computeStressOnQuad
(
grad_u_rot
,
sigma_rot
,
sigma_th
);
auto
sigma
=
this
->
reverseRotation
(
sigma_rot
,
rotation_matrix
);
Matrix
<
Real
>
identity
(
3
,
3
);
identity
.
eye
();
Matrix
<
Real
>
sigma_expected
=
0.5
*
E
/
(
1
+
nu
)
*
(
grad_u
+
grad_u
.
transpose
())
+
sigma_th
*
identity
;
auto
diff
=
sigma
-
sigma_expected
;
Real
stress_error
=
diff
.
norm
<
L_inf
>
();
ASSERT_NEAR
(
stress_error
,
0.
,
1e-14
);
}
/* -------------------------------------------------------------------------- */
template
<>
void
FriendMaterial
<
MaterialElastic
<
3
>>::
testEnergyDensity
()
{
Matrix
<
Real
>
sigma
=
{{
1
,
2
,
3
},
{
2
,
4
,
5
},
{
3
,
5
,
6
}};
Matrix
<
Real
>
eps
=
{{
1
,
0
,
0
},
{
0
,
1
,
0
},
{
0
,
0
,
1
}};
Real
epot
=
0
;
Real
solution
=
5.5
;
this
->
computePotentialEnergyOnQuad
(
eps
,
sigma
,
epot
);
ASSERT_NEAR
(
epot
,
solution
,
1e-14
);
}
/* -------------------------------------------------------------------------- */
template
<>
void
FriendMaterial
<
MaterialElastic
<
3
>>::
testComputeTangentModuli
()
{
Real
E
=
1.
;
Real
nu
=
.3
;
setParam
(
"E"
,
E
);
setParam
(
"nu"
,
nu
);
Matrix
<
Real
>
tangent
(
6
,
6
);
// clang-format off
Matrix
<
Real
>
solution
=
{
{
1
-
nu
,
nu
,
nu
,
0
,
0
,
0
},
{
nu
,
1
-
nu
,
nu
,
0
,
0
,
0
},
{
nu
,
nu
,
1
-
nu
,
0
,
0
,
0
},
{
0
,
0
,
0
,
(
1
-
2
*
nu
)
/
2
,
0
,
0
},
{
0
,
0
,
0
,
0
,
(
1
-
2
*
nu
)
/
2
,
0
},
{
0
,
0
,
0
,
0
,
0
,
(
1
-
2
*
nu
)
/
2
},
};
// clang-format on
solution
*=
E
/
((
1
+
nu
)
*
(
1
-
2
*
nu
));
this
->
computeTangentModuliOnQuad
(
tangent
);
Real
tangent_error
=
(
tangent
-
solution
).
norm
<
L_2
>
();
ASSERT_NEAR
(
tangent_error
,
0
,
1e-14
);
}
/* -------------------------------------------------------------------------- */
template
<>
void
FriendMaterial
<
MaterialElastic
<
3
>>::
testCelerity
()
{
Real
E
=
1.
;
Real
nu
=
.3
;
Real
rho
=
2
;
setParam
(
"E"
,
E
);
setParam
(
"nu"
,
nu
);
setParam
(
"rho"
,
rho
);
auto
push_wave_speed
=
this
->
getPushWaveSpeed
(
Element
());
auto
celerity
=
this
->
getCelerity
(
Element
());
Real
K
=
E
/
(
3
*
(
1
-
2
*
nu
));
Real
mu
=
E
/
(
2
*
(
1
+
nu
));
Real
sol
=
std
::
sqrt
((
K
+
4.
/
3
*
mu
)
/
rho
);
ASSERT_NEAR
(
push_wave_speed
,
sol
,
1e-14
);
ASSERT_NEAR
(
celerity
,
sol
,
1e-14
);
auto
shear_wave_speed
=
this
->
getShearWaveSpeed
(
Element
());
sol
=
std
::
sqrt
(
mu
/
rho
);
ASSERT_NEAR
(
shear_wave_speed
,
sol
,
1e-14
);
}
/* -------------------------------------------------------------------------- */
template
<>
void
FriendMaterial
<
MaterialElasticOrthotropic
<
2
>>::
testComputeStress
()
{
UInt
Dim
=
2
;
Real
E1
=
1.
;
Real
E2
=
2.
;
Real
nu12
=
0.1
;
Real
G12
=
2.
;
setParamNoUpdate
(
"E1"
,
E1
);
setParamNoUpdate
(
"E2"
,
E2
);
setParamNoUpdate
(
"nu12"
,
nu12
);
setParamNoUpdate
(
"G12"
,
G12
);
// material frame of reference is rotate by rotation_matrix starting from
// canonical basis
Matrix
<
Real
>
rotation_matrix
=
getRandomRotation2d
();
// canonical basis as expressed in the material frame of reference, as
// required by MaterialElasticOrthotropic class (it is simply given by the
// columns of the rotation_matrix; the lines give the material basis expressed
// in the canonical frame of reference)
Vector
<
Real
>
v1
(
Dim
);
Vector
<
Real
>
v2
(
Dim
);
for
(
UInt
i
=
0
;
i
<
Dim
;
++
i
)
{
v1
[
i
]
=
rotation_matrix
(
i
,
0
);
v2
[
i
]
=
rotation_matrix
(
i
,
1
);
}
setParamNoUpdate
(
"n1"
,
v1
.
normalize
());
setParamNoUpdate
(
"n2"
,
v2
.
normalize
());
// set internal Cijkl matrix expressed in the canonical frame of reference
this
->
updateInternalParameters
();
// gradient in material frame of reference
auto
grad_u
=
(
this
->
getDeviatoricStrain
(
2.
)
+
this
->
getHydrostaticStrain
(
1.
))
.
block
(
0
,
0
,
2
,
2
);
// gradient in canonical basis (we need to rotate *back* to the canonical
// basis)
auto
grad_u_rot
=
this
->
reverseRotation
(
grad_u
,
rotation_matrix
);
// stress in the canonical basis
Matrix
<
Real
>
sigma_rot
(
2
,
2
);
this
->
computeStressOnQuad
(
grad_u_rot
,
sigma_rot
);
// stress in the material reference (we need to apply the rotation)
auto
sigma
=
this
->
applyRotation
(
sigma_rot
,
rotation_matrix
);
// construction of Cijkl engineering tensor in the *material* frame of
// reference
// ref: http://solidmechanics.org/Text/Chapter3_2/Chapter3_2.php#Sect3_2_13
Real
nu21
=
nu12
*
E2
/
E1
;
Real
gamma
=
1
/
(
1
-
nu12
*
nu21
);
Matrix
<
Real
>
C_expected
(
2
*
Dim
,
2
*
Dim
,
0
);
C_expected
(
0
,
0
)
=
gamma
*
E1
;
C_expected
(
1
,
1
)
=
gamma
*
E2
;
C_expected
(
2
,
2
)
=
G12
;
C_expected
(
1
,
0
)
=
C_expected
(
0
,
1
)
=
gamma
*
E1
*
nu21
;
// epsilon is computed directly in the *material* frame of reference
Matrix
<
Real
>
epsilon
=
0.5
*
(
grad_u
+
grad_u
.
transpose
());
// sigma_expected is computed directly in the *material* frame of reference
Matrix
<
Real
>
sigma_expected
(
Dim
,
Dim
);
for
(
UInt
i
=
0
;
i
<
Dim
;
++
i
)
{
for
(
UInt
j
=
0
;
j
<
Dim
;
++
j
)
{
sigma_expected
(
i
,
i
)
+=
C_expected
(
i
,
j
)
*
epsilon
(
j
,
j
);
}
}
sigma_expected
(
0
,
1
)
=
sigma_expected
(
1
,
0
)
=
C_expected
(
2
,
2
)
*
2
*
epsilon
(
0
,
1
);
// sigmas are checked in the *material* frame of reference
auto
diff
=
sigma
-
sigma_expected
;
Real
stress_error
=
diff
.
norm
<
L_inf
>
();
ASSERT_NEAR
(
stress_error
,
0.
,
1e-13
);
}
/* -------------------------------------------------------------------------- */
template
<>
void
FriendMaterial
<
MaterialElasticOrthotropic
<
2
>>::
testEnergyDensity
()
{
Matrix
<
Real
>
sigma
=
{{
1
,
2
},
{
2
,
4
}};
Matrix
<
Real
>
eps
=
{{
1
,
0
},
{
0
,
1
}};
Real
epot
=
0
;
Real
solution
=
2.5
;
this
->
computePotentialEnergyOnQuad
(
eps
,
sigma
,
epot
);
ASSERT_NEAR
(
epot
,
solution
,
1e-14
);
}
/* -------------------------------------------------------------------------- */
template
<>
void
FriendMaterial
<
MaterialElasticOrthotropic
<
2
>>::
testComputeTangentModuli
()
{
// Note: for this test material and canonical basis coincide
Vector
<
Real
>
n1
=
{
1
,
0
};
Vector
<
Real
>
n2
=
{
0
,
1
};
Real
E1
=
1.
;
Real
E2
=
2.
;
Real
nu12
=
0.1
;
Real
G12
=
2.
;
setParamNoUpdate
(
"n1"
,
n1
);
setParamNoUpdate
(
"n2"
,
n2
);
setParamNoUpdate
(
"E1"
,
E1
);
setParamNoUpdate
(
"E2"
,
E2
);
setParamNoUpdate
(
"nu12"
,
nu12
);
setParamNoUpdate
(
"G12"
,
G12
);
// set internal Cijkl matrix expressed in the canonical frame of reference
this
->
updateInternalParameters
();
// construction of Cijkl engineering tensor in the *material* frame of
// reference
// ref: http://solidmechanics.org/Text/Chapter3_2/Chapter3_2.php#Sect3_2_13
Real
nu21
=
nu12
*
E2
/
E1
;
Real
gamma
=
1
/
(
1
-
nu12
*
nu21
);
Matrix
<
Real
>
C_expected
(
3
,
3
);
C_expected
(
0
,
0
)
=
gamma
*
E1
;
C_expected
(
1
,
1
)
=
gamma
*
E2
;
C_expected
(
2
,
2
)
=
G12
;
C_expected
(
1
,
0
)
=
C_expected
(
0
,
1
)
=
gamma
*
E1
*
nu21
;
Matrix
<
Real
>
tangent
(
3
,
3
);
this
->
computeTangentModuliOnQuad
(
tangent
);
Real
tangent_error
=
(
tangent
-
C_expected
).
norm
<
L_2
>
();
ASSERT_NEAR
(
tangent_error
,
0
,
1e-14
);
}
/* -------------------------------------------------------------------------- */
template
<>
void
FriendMaterial
<
MaterialElasticOrthotropic
<
3
>>::
testComputeStress
()
{
UInt
Dim
=
3
;
Real
E1
=
1.
;
Real
E2
=
2.
;
Real
E3
=
3.
;
Real
nu12
=
0.1
;
Real
nu13
=
0.2
;
Real
nu23
=
0.3
;
Real
G12
=
2.
;
Real
G13
=
3.
;
Real
G23
=
1.
;
setParamNoUpdate
(
"E1"
,
E1
);
setParamNoUpdate
(
"E2"
,
E2
);
setParamNoUpdate
(
"E3"
,
E3
);
setParamNoUpdate
(
"nu12"
,
nu12
);
setParamNoUpdate
(
"nu13"
,
nu13
);
setParamNoUpdate
(
"nu23"
,
nu23
);
setParamNoUpdate
(
"G12"
,
G12
);
setParamNoUpdate
(
"G13"
,
G13
);
setParamNoUpdate
(
"G23"
,
G23
);
// material frame of reference is rotate by rotation_matrix starting from
// canonical basis
Matrix
<
Real
>
rotation_matrix
=
getRandomRotation3d
();
// canonical basis as expressed in the material frame of reference, as
// required by MaterialElasticOrthotropic class (it is simply given by the
// columns of the rotation_matrix; the lines give the material basis expressed
// in the canonical frame of reference)
Vector
<
Real
>
v1
(
Dim
);
Vector
<
Real
>
v2
(
Dim
);
Vector
<
Real
>
v3
(
Dim
);
for
(
UInt
i
=
0
;
i
<
Dim
;
++
i
)
{
v1
[
i
]
=
rotation_matrix
(
i
,
0
);
v2
[
i
]
=
rotation_matrix
(
i
,
1
);
v3
[
i
]
=
rotation_matrix
(
i
,
2
);
}
setParamNoUpdate
(
"n1"
,
v1
.
normalize
());
setParamNoUpdate
(
"n2"
,
v2
.
normalize
());
setParamNoUpdate
(
"n3"
,
v3
.
normalize
());
// set internal Cijkl matrix expressed in the canonical frame of reference
this
->
updateInternalParameters
();
// gradient in material frame of reference
auto
grad_u
=
this
->
getDeviatoricStrain
(
2.
)
+
this
->
getHydrostaticStrain
(
1.
);
// gradient in canonical basis (we need to rotate *back* to the canonical
// basis)
auto
grad_u_rot
=
this
->
reverseRotation
(
grad_u
,
rotation_matrix
);
// stress in the canonical basis
Matrix
<
Real
>
sigma_rot
(
3
,
3
);
this
->
computeStressOnQuad
(
grad_u_rot
,
sigma_rot
);
// stress in the material reference (we need to apply the rotation)
auto
sigma
=
this
->
applyRotation
(
sigma_rot
,
rotation_matrix
);
// construction of Cijkl engineering tensor in the *material* frame of
// reference
// ref: http://solidmechanics.org/Text/Chapter3_2/Chapter3_2.php#Sect3_2_13
Real
nu21
=
nu12
*
E2
/
E1
;
Real
nu31
=
nu13
*
E3
/
E1
;
Real
nu32
=
nu23
*
E3
/
E2
;
Real
gamma
=
1
/
(
1
-
nu12
*
nu21
-
nu23
*
nu32
-
nu31
*
nu13
-
2
*
nu21
*
nu32
*
nu13
);
Matrix
<
Real
>
C_expected
(
6
,
6
);
C_expected
(
0
,
0
)
=
gamma
*
E1
*
(
1
-
nu23
*
nu32
);
C_expected
(
1
,
1
)
=
gamma
*
E2
*
(
1
-
nu13
*
nu31
);
C_expected
(
2
,
2
)
=
gamma
*
E3
*
(
1
-
nu12
*
nu21
);
C_expected
(
1
,
0
)
=
C_expected
(
0
,
1
)
=
gamma
*
E1
*
(
nu21
+
nu31
*
nu23
);
C_expected
(
2
,
0
)
=
C_expected
(
0
,
2
)
=
gamma
*
E1
*
(
nu31
+
nu21
*
nu32
);
C_expected
(
2
,
1
)
=
C_expected
(
1
,
2
)
=
gamma
*
E2
*
(
nu32
+
nu12
*
nu31
);
C_expected
(
3
,
3
)
=
G23
;
C_expected
(
4
,
4
)
=
G13
;
C_expected
(
5
,
5
)
=
G12
;
// epsilon is computed directly in the *material* frame of reference
Matrix
<
Real
>
epsilon
=
0.5
*
(
grad_u
+
grad_u
.
transpose
());
// sigma_expected is computed directly in the *material* frame of reference
Matrix
<
Real
>
sigma_expected
(
Dim
,
Dim
);
for
(
UInt
i
=
0
;
i
<
Dim
;
++
i
)
{
for
(
UInt
j
=
0
;
j
<
Dim
;
++
j
)
{
sigma_expected
(
i
,
i
)
+=
C_expected
(
i
,
j
)
*
epsilon
(
j
,
j
);
}
}
sigma_expected
(
0
,
1
)
=
C_expected
(
5
,
5
)
*
2
*
epsilon
(
0
,
1
);
sigma_expected
(
0
,
2
)
=
C_expected
(
4
,
4
)
*
2
*
epsilon
(
0
,
2
);
sigma_expected
(
1
,
2
)
=
C_expected
(
3
,
3
)
*
2
*
epsilon
(
1
,
2
);
sigma_expected
(
1
,
0
)
=
sigma_expected
(
0
,
1
);
sigma_expected
(
2
,
0
)
=
sigma_expected
(
0
,
2
);
sigma_expected
(
2
,
1
)
=
sigma_expected
(
1
,
2
);
// sigmas are checked in the *material* frame of reference
auto
diff
=
sigma
-
sigma_expected
;
Real
stress_error
=
diff
.
norm
<
L_inf
>
();
ASSERT_NEAR
(
stress_error
,
0.
,
1e-13
);
}
/* -------------------------------------------------------------------------- */
template
<>
void
FriendMaterial
<
MaterialElasticOrthotropic
<
3
>>::
testEnergyDensity
()
{
Matrix
<
Real
>
sigma
=
{{
1
,
2
,
3
},
{
2
,
4
,
5
},
{
3
,
5
,
6
}};
Matrix
<
Real
>
eps
=
{{
1
,
0
,
0
},
{
0
,
1
,
0
},
{
0
,
0
,
1
}};
Real
epot
=
0
;
Real
solution
=
5.5
;
this
->
computePotentialEnergyOnQuad
(
eps
,
sigma
,
epot
);
ASSERT_NEAR
(
epot
,
solution
,
1e-14
);
}
/* -------------------------------------------------------------------------- */
template
<>
void
FriendMaterial
<
MaterialElasticOrthotropic
<
3
>>::
testComputeTangentModuli
()
{
// Note: for this test material and canonical basis coincide
UInt
Dim
=
3
;
Vector
<
Real
>
n1
=
{
1
,
0
,
0
};
Vector
<
Real
>
n2
=
{
0
,
1
,
0
};
Vector
<
Real
>
n3
=
{
0
,
0
,
1
};
Real
E1
=
1.
;
Real
E2
=
2.
;
Real
E3
=
3.
;
Real
nu12
=
0.1
;
Real
nu13
=
0.2
;
Real
nu23
=
0.3
;
Real
G12
=
2.
;
Real
G13
=
3.
;
Real
G23
=
1.
;
setParamNoUpdate
(
"n1"
,
n1
);
setParamNoUpdate
(
"n2"
,
n2
);
setParamNoUpdate
(
"n3"
,
n3
);
setParamNoUpdate
(
"E1"
,
E1
);
setParamNoUpdate
(
"E2"
,
E2
);
setParamNoUpdate
(
"E3"
,
E3
);
setParamNoUpdate
(
"nu12"
,
nu12
);
setParamNoUpdate
(
"nu13"
,
nu13
);
setParamNoUpdate
(
"nu23"
,
nu23
);
setParamNoUpdate
(
"G12"
,
G12
);
setParamNoUpdate
(
"G13"
,
G13
);
setParamNoUpdate
(
"G23"
,
G23
);
// set internal Cijkl matrix expressed in the canonical frame of reference
this
->
updateInternalParameters
();
// construction of Cijkl engineering tensor in the *material* frame of
// reference
// ref: http://solidmechanics.org/Text/Chapter3_2/Chapter3_2.php#Sect3_2_13
Real
nu21
=
nu12
*
E2
/
E1
;
Real
nu31
=
nu13
*
E3
/
E1
;
Real
nu32
=
nu23
*
E3
/
E2
;
Real
gamma
=
1
/
(
1
-
nu12
*
nu21
-
nu23
*
nu32
-
nu31
*
nu13
-
2
*
nu21
*
nu32
*
nu13
);
Matrix
<
Real
>
C_expected
(
2
*
Dim
,
2
*
Dim
,
0
);
C_expected
(
0
,
0
)
=
gamma
*
E1
*
(
1
-
nu23
*
nu32
);
C_expected
(
1
,
1
)
=
gamma
*
E2
*
(
1
-
nu13
*
nu31
);
C_expected
(
2
,
2
)
=
gamma
*
E3
*
(
1
-
nu12
*
nu21
);
C_expected
(
1
,
0
)
=
C_expected
(
0
,
1
)
=
gamma
*
E1
*
(
nu21
+
nu31
*
nu23
);
C_expected
(
2
,
0
)
=
C_expected
(
0
,
2
)
=
gamma
*
E1
*
(
nu31
+
nu21
*
nu32
);
C_expected
(
2
,
1
)
=
C_expected
(
1
,
2
)
=
gamma
*
E2
*
(
nu32
+
nu12
*
nu31
);
C_expected
(
3
,
3
)
=
G23
;
C_expected
(
4
,
4
)
=
G13
;
C_expected
(
5
,
5
)
=
G12
;
Matrix
<
Real
>
tangent
(
6
,
6
);
this
->
computeTangentModuliOnQuad
(
tangent
);
Real
tangent_error
=
(
tangent
-
C_expected
).
norm
<
L_2
>
();
ASSERT_NEAR
(
tangent_error
,
0
,
1e-14
);
}
/* -------------------------------------------------------------------------- */
template
<>
void
FriendMaterial
<
MaterialElasticLinearAnisotropic
<
2
>>::
testComputeStress
()
{
UInt
Dim
=
2
;
Matrix
<
Real
>
C
=
{
{
1.0
,
0.3
,
0.4
},
{
0.3
,
2.0
,
0.1
},
{
0.4
,
0.1
,
1.5
},
};
setParamNoUpdate
(
"C11"
,
C
(
0
,
0
));
setParamNoUpdate
(
"C12"
,
C
(
0
,
1
));
setParamNoUpdate
(
"C13"
,
C
(
0
,
2
));
setParamNoUpdate
(
"C22"
,
C
(
1
,
1
));
setParamNoUpdate
(
"C23"
,
C
(
1
,
2
));
setParamNoUpdate
(
"C33"
,
C
(
2
,
2
));
// material frame of reference is rotate by rotation_matrix starting from
// canonical basis
Matrix
<
Real
>
rotation_matrix
=
getRandomRotation2d
();
// canonical basis as expressed in the material frame of reference, as
// required by MaterialElasticLinearAnisotropic class (it is simply given by
// the columns of the rotation_matrix; the lines give the material basis
// expressed in the canonical frame of reference)
Vector
<
Real
>
v1
(
Dim
);
Vector
<
Real
>
v2
(
Dim
);
for
(
UInt
i
=
0
;
i
<
Dim
;
++
i
)
{
v1
[
i
]
=
rotation_matrix
(
i
,
0
);
v2
[
i
]
=
rotation_matrix
(
i
,
1
);
}
setParamNoUpdate
(
"n1"
,
v1
.
normalize
());
setParamNoUpdate
(
"n2"
,
v2
.
normalize
());
// set internal Cijkl matrix expressed in the canonical frame of reference
this
->
updateInternalParameters
();
// gradient in material frame of reference
auto
grad_u
=
(
this
->
getDeviatoricStrain
(
2.
)
+
this
->
getHydrostaticStrain
(
1.
))
.
block
(
0
,
0
,
2
,
2
);
// gradient in canonical basis (we need to rotate *back* to the canonical
// basis)
auto
grad_u_rot
=
this
->
reverseRotation
(
grad_u
,
rotation_matrix
);
// stress in the canonical basis
Matrix
<
Real
>
sigma_rot
(
2
,
2
);
this
->
computeStressOnQuad
(
grad_u_rot
,
sigma_rot
);
// stress in the material reference (we need to apply the rotation)
auto
sigma
=
this
->
applyRotation
(
sigma_rot
,
rotation_matrix
);
// epsilon is computed directly in the *material* frame of reference
Matrix
<
Real
>
epsilon
=
0.5
*
(
grad_u
+
grad_u
.
transpose
());
Vector
<
Real
>
epsilon_voigt
(
3
);
epsilon_voigt
(
0
)
=
epsilon
(
0
,
0
);
epsilon_voigt
(
1
)
=
epsilon
(
1
,
1
);
epsilon_voigt
(
2
)
=
2
*
epsilon
(
0
,
1
);
// sigma_expected is computed directly in the *material* frame of reference
Vector
<
Real
>
sigma_voigt
=
C
*
epsilon_voigt
;
Matrix
<
Real
>
sigma_expected
(
Dim
,
Dim
);
sigma_expected
(
0
,
0
)
=
sigma_voigt
(
0
);
sigma_expected
(
1
,
1
)
=
sigma_voigt
(
1
);
sigma_expected
(
0
,
1
)
=
sigma_expected
(
1
,
0
)
=
sigma_voigt
(
2
);
// sigmas are checked in the *material* frame of reference
auto
diff
=
sigma
-
sigma_expected
;
Real
stress_error
=
diff
.
norm
<
L_inf
>
();
ASSERT_NEAR
(
stress_error
,
0.
,
1e-13
);
}
/* -------------------------------------------------------------------------- */
template
<>
void
FriendMaterial
<
MaterialElasticLinearAnisotropic
<
2
>>::
testEnergyDensity
()
{
Matrix
<
Real
>
sigma
=
{{
1
,
2
},
{
2
,
4
}};
Matrix
<
Real
>
eps
=
{{
1
,
0
},
{
0
,
1
}};
Real
epot
=
0
;
Real
solution
=
2.5
;
this
->
computePotentialEnergyOnQuad
(
eps
,
sigma
,
epot
);
ASSERT_NEAR
(
epot
,
solution
,
1e-14
);
}
/* -------------------------------------------------------------------------- */
template
<>
void
FriendMaterial
<
MaterialElasticLinearAnisotropic
<
2
>>::
testComputeTangentModuli
()
{
// Note: for this test material and canonical basis coincide
Matrix
<
Real
>
C
=
{
{
1.0
,
0.3
,
0.4
},
{
0.3
,
2.0
,
0.1
},
{
0.4
,
0.1
,
1.5
},
};
setParamNoUpdate
(
"C11"
,
C
(
0
,
0
));
setParamNoUpdate
(
"C12"
,
C
(
0
,
1
));
setParamNoUpdate
(
"C13"
,
C
(
0
,
2
));
setParamNoUpdate
(
"C22"
,
C
(
1
,
1
));
setParamNoUpdate
(
"C23"
,
C
(
1
,
2
));
setParamNoUpdate
(
"C33"
,
C
(
2
,
2
));
// set internal Cijkl matrix expressed in the canonical frame of reference
this
->
updateInternalParameters
();
Matrix
<
Real
>
tangent
(
3
,
3
);
this
->
computeTangentModuliOnQuad
(
tangent
);
Real
tangent_error
=
(
tangent
-
C
).
norm
<
L_2
>
();
ASSERT_NEAR
(
tangent_error
,
0
,
1e-14
);
}
/* -------------------------------------------------------------------------- */
template
<>
void
FriendMaterial
<
MaterialElasticLinearAnisotropic
<
3
>>::
testComputeStress
()
{
UInt
Dim
=
3
;
Matrix
<
Real
>
C
=
{
{
1.0
,
0.3
,
0.4
,
0.3
,
0.2
,
0.1
},
{
0.3
,
2.0
,
0.1
,
0.2
,
0.3
,
0.2
},
{
0.4
,
0.1
,
1.5
,
0.1
,
0.4
,
0.3
},
{
0.3
,
0.2
,
0.1
,
2.4
,
0.1
,
0.4
},
{
0.2
,
0.3
,
0.4
,
0.1
,
0.9
,
0.1
},
{
0.1
,
0.2
,
0.3
,
0.4
,
0.1
,
1.2
},
};
setParamNoUpdate
(
"C11"
,
C
(
0
,
0
));
setParamNoUpdate
(
"C12"
,
C
(
0
,
1
));
setParamNoUpdate
(
"C13"
,
C
(
0
,
2
));
setParamNoUpdate
(
"C14"
,
C
(
0
,
3
));
setParamNoUpdate
(
"C15"
,
C
(
0
,
4
));
setParamNoUpdate
(
"C16"
,
C
(
0
,
5
));
setParamNoUpdate
(
"C22"
,
C
(
1
,
1
));
setParamNoUpdate
(
"C23"
,
C
(
1
,
2
));
setParamNoUpdate
(
"C24"
,
C
(
1
,
3
));
setParamNoUpdate
(
"C25"
,
C
(
1
,
4
));
setParamNoUpdate
(
"C26"
,
C
(
1
,
5
));
setParamNoUpdate
(
"C33"
,
C
(
2
,
2
));
setParamNoUpdate
(
"C34"
,
C
(
2
,
3
));
setParamNoUpdate
(
"C35"
,
C
(
2
,
4
));
setParamNoUpdate
(
"C36"
,
C
(
2
,
5
));
setParamNoUpdate
(
"C44"
,
C
(
3
,
3
));
setParamNoUpdate
(
"C45"
,
C
(
3
,
4
));
setParamNoUpdate
(
"C46"
,
C
(
3
,
5
));
setParamNoUpdate
(
"C55"
,
C
(
4
,
4
));
setParamNoUpdate
(
"C56"
,
C
(
4
,
5
));
setParamNoUpdate
(
"C66"
,
C
(
5
,
5
));
// material frame of reference is rotate by rotation_matrix starting from
// canonical basis
Matrix
<
Real
>
rotation_matrix
=
getRandomRotation3d
();
// canonical basis as expressed in the material frame of reference, as
// required by MaterialElasticLinearAnisotropic class (it is simply given by
// the columns of the rotation_matrix; the lines give the material basis
// expressed in the canonical frame of reference)
Vector
<
Real
>
v1
(
Dim
);
Vector
<
Real
>
v2
(
Dim
);
Vector
<
Real
>
v3
(
Dim
);
for
(
UInt
i
=
0
;
i
<
Dim
;
++
i
)
{
v1
[
i
]
=
rotation_matrix
(
i
,
0
);
v2
[
i
]
=
rotation_matrix
(
i
,
1
);
v3
[
i
]
=
rotation_matrix
(
i
,
2
);
}
setParamNoUpdate
(
"n1"
,
v1
.
normalize
());
setParamNoUpdate
(
"n2"
,
v2
.
normalize
());
setParamNoUpdate
(
"n3"
,
v3
.
normalize
());
// set internal Cijkl matrix expressed in the canonical frame of reference
this
->
updateInternalParameters
();
// gradient in material frame of reference
auto
grad_u
=
this
->
getDeviatoricStrain
(
2.
)
+
this
->
getHydrostaticStrain
(
1.
);
// gradient in canonical basis (we need to rotate *back* to the canonical
// basis)
auto
grad_u_rot
=
this
->
reverseRotation
(
grad_u
,
rotation_matrix
);
// stress in the canonical basis
Matrix
<
Real
>
sigma_rot
(
3
,
3
);
this
->
computeStressOnQuad
(
grad_u_rot
,
sigma_rot
);
// stress in the material reference (we need to apply the rotation)
auto
sigma
=
this
->
applyRotation
(
sigma_rot
,
rotation_matrix
);
// epsilon is computed directly in the *material* frame of reference
Matrix
<
Real
>
epsilon
=
0.5
*
(
grad_u
+
grad_u
.
transpose
());
Vector
<
Real
>
epsilon_voigt
(
6
);
epsilon_voigt
(
0
)
=
epsilon
(
0
,
0
);
epsilon_voigt
(
1
)
=
epsilon
(
1
,
1
);
epsilon_voigt
(
2
)
=
epsilon
(
2
,
2
);
epsilon_voigt
(
3
)
=
2
*
epsilon
(
1
,
2
);
epsilon_voigt
(
4
)
=
2
*
epsilon
(
0
,
2
);
epsilon_voigt
(
5
)
=
2
*
epsilon
(
0
,
1
);
// sigma_expected is computed directly in the *material* frame of reference
Vector
<
Real
>
sigma_voigt
=
C
*
epsilon_voigt
;
Matrix
<
Real
>
sigma_expected
(
Dim
,
Dim
);
sigma_expected
(
0
,
0
)
=
sigma_voigt
(
0
);
sigma_expected
(
1
,
1
)
=
sigma_voigt
(
1
);
sigma_expected
(
2
,
2
)
=
sigma_voigt
(
2
);
sigma_expected
(
1
,
2
)
=
sigma_expected
(
2
,
1
)
=
sigma_voigt
(
3
);
sigma_expected
(
0
,
2
)
=
sigma_expected
(
2
,
0
)
=
sigma_voigt
(
4
);
sigma_expected
(
0
,
1
)
=
sigma_expected
(
1
,
0
)
=
sigma_voigt
(
5
);
// sigmas are checked in the *material* frame of reference
auto
diff
=
sigma
-
sigma_expected
;
Real
stress_error
=
diff
.
norm
<
L_inf
>
();
ASSERT_NEAR
(
stress_error
,
0.
,
1e-13
);
}
/* -------------------------------------------------------------------------- */
template
<>
void
FriendMaterial
<
MaterialElasticLinearAnisotropic
<
3
>>::
testEnergyDensity
()
{
Matrix
<
Real
>
sigma
=
{{
1
,
2
,
3
},
{
2
,
4
,
5
},
{
3
,
5
,
6
}};
Matrix
<
Real
>
eps
=
{{
1
,
0
,
0
},
{
0
,
1
,
0
},
{
0
,
0
,
1
}};
Real
epot
=
0
;
Real
solution
=
5.5
;
this
->
computePotentialEnergyOnQuad
(
eps
,
sigma
,
epot
);
ASSERT_NEAR
(
epot
,
solution
,
1e-14
);
}
/* -------------------------------------------------------------------------- */
template
<>
void
FriendMaterial
<
MaterialElasticLinearAnisotropic
<
3
>>::
testComputeTangentModuli
()
{
// Note: for this test material and canonical basis coincide
Matrix
<
Real
>
C
=
{
{
1.0
,
0.3
,
0.4
,
0.3
,
0.2
,
0.1
},
{
0.3
,
2.0
,
0.1
,
0.2
,
0.3
,
0.2
},
{
0.4
,
0.1
,
1.5
,
0.1
,
0.4
,
0.3
},
{
0.3
,
0.2
,
0.1
,
2.4
,
0.1
,
0.4
},
{
0.2
,
0.3
,
0.4
,
0.1
,
0.9
,
0.1
},
{
0.1
,
0.2
,
0.3
,
0.4
,
0.1
,
1.2
},
};
setParamNoUpdate
(
"C11"
,
C
(
0
,
0
));
setParamNoUpdate
(
"C12"
,
C
(
0
,
1
));
setParamNoUpdate
(
"C13"
,
C
(
0
,
2
));
setParamNoUpdate
(
"C14"
,
C
(
0
,
3
));
setParamNoUpdate
(
"C15"
,
C
(
0
,
4
));
setParamNoUpdate
(
"C16"
,
C
(
0
,
5
));
setParamNoUpdate
(
"C22"
,
C
(
1
,
1
));
setParamNoUpdate
(
"C23"
,
C
(
1
,
2
));
setParamNoUpdate
(
"C24"
,
C
(
1
,
3
));
setParamNoUpdate
(
"C25"
,
C
(
1
,
4
));
setParamNoUpdate
(
"C26"
,
C
(
1
,
5
));
setParamNoUpdate
(
"C33"
,
C
(
2
,
2
));
setParamNoUpdate
(
"C34"
,
C
(
2
,
3
));
setParamNoUpdate
(
"C35"
,
C
(
2
,
4
));
setParamNoUpdate
(
"C36"
,
C
(
2
,
5
));
setParamNoUpdate
(
"C44"
,
C
(
3
,
3
));
setParamNoUpdate
(
"C45"
,
C
(
3
,
4
));
setParamNoUpdate
(
"C46"
,
C
(
3
,
5
));
setParamNoUpdate
(
"C55"
,
C
(
4
,
4
));
setParamNoUpdate
(
"C56"
,
C
(
4
,
5
));
setParamNoUpdate
(
"C66"
,
C
(
5
,
5
));
// set internal Cijkl matrix expressed in the canonical frame of reference
this
->
updateInternalParameters
();
Matrix
<
Real
>
tangent
(
6
,
6
);
this
->
computeTangentModuliOnQuad
(
tangent
);
Real
tangent_error
=
(
tangent
-
C
).
norm
<
L_2
>
();
ASSERT_NEAR
(
tangent_error
,
0
,
1e-14
);
}
/* -------------------------------------------------------------------------- */
namespace
{
template
<
typename
T
>
class
TestElasticMaterialFixture
:
public
::
TestMaterialFixture
<
T
>
{};
TYPED_TEST_CASE
(
TestElasticMaterialFixture
,
types
);
TYPED_TEST
(
TestElasticMaterialFixture
,
ElasticComputeStress
)
{
this
->
material
->
testComputeStress
();
}
TYPED_TEST
(
TestElasticMaterialFixture
,
ElasticEnergyDensity
)
{
this
->
material
->
testEnergyDensity
();
}
TYPED_TEST
(
TestElasticMaterialFixture
,
ElasticComputeTangentModuli
)
{
this
->
material
->
testComputeTangentModuli
();
}
TYPED_TEST
(
TestElasticMaterialFixture
,
ElasticComputeCelerity
)
{
this
->
material
->
testCelerity
();
}
}
// namespace
Event Timeline
Log In to Comment