Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F78846470
test_integral_operators.py
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, Aug 23, 07:40
Size
7 KB
Mime Type
text/x-python
Expires
Sun, Aug 25, 07:40 (1 d, 22 h)
Engine
blob
Format
Raw Data
Handle
20102399
Attached To
rTAMAAS tamaas
test_integral_operators.py
View Options
# -*- coding: utf-8 -*-
# @file
# @section LICENSE
#
# Copyright (©) 2016-2021 EPFL (École Polytechnique Fédérale de Lausanne),
# Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU Affero General Public License as published
# by the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program 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 Affero General Public License for more details.
#
# You should have received a copy of the GNU Affero General Public License
# along with this program. If not, see <https://www.gnu.org/licenses/>.
import
tamaas
as
tm
import
numpy
as
np
import
pytest
from
numpy.linalg
import
norm
from
register_integral_operators
import
register_kelvin_force
,
\
set_integration_method
@pytest.fixture
(
scope
=
"module"
,
params
=
[
tm
.
integration_method
.
linear
,
tm
.
integration_method
.
cutoff
])
def
integral_fixture
(
tamaas_fixture
,
request
):
return
request
.
param
def
test_kelvin_volume_force
(
integral_fixture
):
N
=
65
E
=
1.
nu
=
0.3
mu
=
E
/
(
2
*
(
1
+
nu
))
domain
=
np
.
array
([
1.
]
*
3
)
omega
=
2
*
np
.
pi
*
np
.
array
([
1
,
1
])
/
domain
[:
2
]
omega_
=
norm
(
omega
)
discretization
=
[
N
]
*
3
model
=
tm
.
ModelFactory
.
createModel
(
tm
.
model_type
.
volume_2d
,
domain
,
discretization
)
model
.
E
=
E
model
.
nu
=
nu
register_kelvin_force
(
model
)
kelvin
=
model
.
operators
[
'kelvin_force'
]
set_integration_method
(
kelvin
,
integral_fixture
,
1e-12
)
coords
=
[
np
.
linspace
(
0
,
domain
[
i
],
discretization
[
i
],
endpoint
=
False
,
dtype
=
tm
.
dtype
)
for
i
in
range
(
2
)]
\
+
[
np
.
linspace
(
0
,
domain
[
2
],
discretization
[
2
])]
x
,
y
=
np
.
meshgrid
(
*
coords
[:
2
],
indexing
=
'ij'
)
displacement
=
model
[
'displacement'
]
source
=
np
.
zeros_like
(
displacement
)
# The integral of forces should stay constant
source
[
N
//
2
,
:,
:,
2
]
=
np
.
sin
(
omega
[
0
]
*
x
)
*
np
.
sin
(
omega
[
1
]
*
y
)
*
(
N
-
1
)
kelvin
(
source
,
displacement
)
z
=
coords
[
2
]
-
0.5
z
,
x
,
y
=
np
.
meshgrid
(
z
,
*
coords
[:
2
],
indexing
=
'ij'
)
solution
=
np
.
zeros_like
(
source
)
solution
[:,
:,
:,
0
]
=
-
np
.
exp
(
-
omega_
*
np
.
abs
(
z
))
/
(
8
*
mu
*
(
1
-
nu
)
*
omega_
)
\
*
omega
[
0
]
*
z
*
np
.
cos
(
omega
[
0
]
*
x
)
*
np
.
sin
(
omega
[
1
]
*
y
)
solution
[:,
:,
:,
1
]
=
-
np
.
exp
(
-
omega_
*
np
.
abs
(
z
))
/
(
8
*
mu
*
(
1
-
nu
)
*
omega_
)
\
*
omega
[
1
]
*
z
*
np
.
sin
(
omega
[
0
]
*
x
)
*
np
.
cos
(
omega
[
1
]
*
y
)
solution
[:,
:,
:,
2
]
=
np
.
exp
(
-
omega_
*
np
.
abs
(
z
))
/
(
8
*
mu
*
(
1
-
nu
)
*
omega_
)
\
*
(
3
-
4
*
nu
+
omega_
*
np
.
abs
(
z
))
*
np
.
sin
(
omega
[
0
]
*
x
)
*
np
.
sin
(
omega
[
1
]
*
y
)
error
=
norm
(
displacement
-
solution
)
/
norm
(
solution
)
assert
error
<
5e-2
def
test_mindlin
(
integral_fixture
):
# Definition of modeled domain
# tm.set_log_level(tm.LogLevel.debug)
model_type
=
tm
.
model_type
.
volume_2d
discretization
=
[
126
,
128
,
128
]
system_size
=
[
1.
,
3.
,
3.
]
integration_method
=
integral_fixture
print
(
integration_method
)
# Material contants
E
=
1.
# Young's modulus
nu
=
0.3
# Poisson's ratio
h
=
0.01
*
E
# Hardening modulus
sigma_y
=
0.3
*
E
# Initial yield stress
# Creation of model
model
=
tm
.
ModelFactory
.
createModel
(
model_type
,
system_size
,
discretization
)
model
.
E
=
E
model
.
nu
=
nu
mu
=
E
/
(
2
*
(
1
+
nu
))
lamda
=
E
*
nu
/
((
1
+
nu
)
*
(
1
-
2
*
nu
))
# Setup for integral operators
residual
=
tm
.
ModelFactory
.
createResidual
(
model
,
sigma_y
,
h
)
residual
.
setIntegrationMethod
(
integration_method
,
1e-12
)
# Coordinates
x
=
np
.
linspace
(
0
,
system_size
[
1
],
discretization
[
1
],
endpoint
=
False
,
dtype
=
tm
.
dtype
)
y
=
np
.
linspace
(
0
,
system_size
[
2
],
discretization
[
2
],
endpoint
=
False
,
dtype
=
tm
.
dtype
)
z
=
np
.
linspace
(
0
,
system_size
[
0
],
discretization
[
0
],
endpoint
=
True
,
dtype
=
tm
.
dtype
)
z
,
x
,
y
=
np
.
meshgrid
(
z
,
x
,
y
,
indexing
=
'ij'
)
# Inclusion definition
a
,
c
=
0.1
,
0.2
center
=
[
system_size
[
1
]
/
2
,
system_size
[
2
]
/
2
,
c
]
r
=
np
.
sqrt
((
x
-
center
[
0
])
**
2
+
(
y
-
center
[
1
])
**
2
+
(
z
-
center
[
2
])
**
2
)
ball
=
r
<
a
# Eigenstrain definition
alpha
=
1
# constant isotropic strain
beta
=
(
3
*
lamda
+
2
*
mu
)
*
alpha
*
np
.
eye
(
3
)
eigenstress
=
np
.
zeros
(
discretization
+
[
3
,
3
],
dtype
=
tm
.
dtype
)
eigenstress
[
ball
,
...
]
=
beta
eigenstress
=
tm
.
compute
.
to_voigt
(
eigenstress
.
reshape
(
discretization
+
[
9
]))
# Array references
gradient
=
residual
.
getVector
()
stress
=
residual
.
getStress
()
# Applying operator
# mindlin = model.operators["mindlin"]
mindlin_gradient
=
model
.
operators
[
"mindlin_gradient"
]
# Not testing displacements yet
# mindlin(eigenstress, model['displacement'])
# Applying gradient
mindlin_gradient
(
eigenstress
,
gradient
)
model
.
operators
[
'hooke'
](
gradient
,
stress
)
stress
-=
eigenstress
# Normalizing stess as in Mura (1976)
T
,
alpha
=
1.
,
1.
beta
=
alpha
*
T
*
(
1
+
nu
)
/
(
1
-
nu
)
stress
*=
1.
/
(
2
*
mu
*
beta
)
n
=
discretization
[
1
]
//
2
vertical_stress
=
stress
[:,
n
,
n
,
:]
z_all
=
np
.
linspace
(
0
,
1
,
discretization
[
0
],
dtype
=
tm
.
dtype
)
sigma_z
=
np
.
zeros_like
(
z_all
)
sigma_t
=
np
.
zeros_like
(
z_all
)
inclusion
=
np
.
abs
(
z_all
-
c
)
<
a
# Computing stresses for exterior points
z
=
z_all
[
~
inclusion
]
R_1
=
np
.
abs
(
z
-
c
)
R_2
=
np
.
abs
(
z
+
c
)
sigma_z
[
~
inclusion
]
=
2
*
mu
*
beta
*
a
**
3
/
3
*
(
1
/
R_1
**
3
-
1
/
R_2
**
3
-
18
*
z
*
(
z
+
c
)
/
R_2
**
5
+
3
*
(
z
+
c
)
**
2
/
R_2
**
5
-
3
*
(
z
-
c
)
**
2
/
R_1
**
5
+
30
*
z
*
(
z
+
c
)
**
3
/
R_2
**
7
)
sigma_t
[
~
inclusion
]
=
2
*
mu
*
beta
*
a
**
3
/
3
*
(
1
/
R_1
**
3
+
(
3
-
8
*
nu
)
/
R_2
**
3
-
6
*
z
*
(
z
+
c
)
/
R_2
**
5
+
12
*
nu
*
(
z
+
c
)
**
2
/
R_2
**
5
)
# Computing stresses for interior points
z
=
z_all
[
inclusion
]
R_1
=
np
.
abs
(
z
-
c
)
R_2
=
np
.
abs
(
z
+
c
)
sigma_z
[
inclusion
]
=
2
*
mu
*
beta
*
a
**
3
/
3
*
(
-
2
/
a
**
3
-
1
/
R_2
**
3
-
18
*
z
*
(
z
+
c
)
/
R_2
**
5
+
3
*
(
z
+
c
)
**
2
/
R_2
**
5
+
30
*
z
*
(
z
+
c
)
**
3
/
R_2
**
7
)
sigma_t
[
inclusion
]
=
2
*
mu
*
beta
*
a
**
3
/
3
*
(
-
2
/
a
**
3
+
(
3
-
8
*
nu
)
/
R_2
**
3
-
6
*
z
*
(
z
+
c
)
/
R_2
**
5
+
12
*
nu
*
(
z
+
c
)
**
2
/
R_2
**
5
)
# This test can be used to debug
if
False
:
import
matplotlib.pyplot
as
plt
plt
.
plot
(
z_all
,
vertical_stress
[:,
2
])
plt
.
plot
(
z_all
,
sigma_z
/
(
2
*
mu
*
beta
))
plt
.
figure
()
plt
.
plot
(
z_all
,
vertical_stress
[:,
0
])
plt
.
plot
(
z_all
,
sigma_t
/
(
2
*
mu
*
beta
))
plt
.
show
()
z_error
=
norm
(
vertical_stress
[:,
2
]
-
sigma_z
/
(
2
*
mu
*
beta
))
\
/
discretization
[
0
]
t_error
=
norm
(
vertical_stress
[:,
0
]
-
sigma_t
/
(
2
*
mu
*
beta
))
\
/
discretization
[
0
]
assert
z_error
<
1e-3
and
t_error
<
1e-3
Event Timeline
Log In to Comment