Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F78665460
bi-material.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
Thu, Aug 22, 07:00
Size
5 KB
Mime Type
text/x-python
Expires
Sat, Aug 24, 07:00 (1 d, 21 h)
Engine
blob
Format
Raw Data
Handle
20050571
Attached To
rAKA akantu
bi-material.py
View Options
__copyright__
=
(
"Copyright (©) 2018-2023 EPFL (Ecole Polytechnique Fédérale de Lausanne)"
"Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)"
)
__license__
=
"LGPLv3"
import
akantu
as
aka
import
numpy
as
np
# ------------------------------------------------------------------------------
class
LocalElastic
(
aka
.
Material
):
def
__init__
(
self
,
model
,
_id
):
super
()
.
__init__
(
model
,
_id
)
super
()
.
registerParamReal
(
'E'
,
aka
.
_pat_readable
|
aka
.
_pat_parsable
,
'Youngs modulus'
)
super
()
.
registerParamReal
(
'nu'
,
aka
.
_pat_readable
|
aka
.
_pat_parsable
,
'Poisson ratio'
)
# change it to have the initialize wrapped
self
.
factor
=
super
()
.
registerInternalReal
(
'factor'
,
1
)
self
.
quad_coords
=
super
()
.
registerInternalReal
(
'quad_coordinates'
,
2
)
def
initMaterial
(
self
):
nu
=
self
.
getReal
(
'nu'
)
E
=
self
.
getReal
(
'E'
)
self
.
mu
=
E
/
(
2
*
(
1
+
nu
))
self
.
lame_lambda
=
nu
*
E
/
(
(
1.
+
nu
)
*
(
1.
-
2.
*
nu
))
# Second Lame coefficient (shear modulus)
self
.
lame_mu
=
E
/
(
2.
*
(
1.
+
nu
))
super
()
.
initMaterial
()
model
=
self
.
getModel
()
model
.
getFEEngine
()
.
computeIntegrationPointsCoordinates
(
self
.
quad_coords
,
self
.
getElementFilter
())
for
elem_type
in
self
.
factor
.
elementTypes
():
factor
=
self
.
factor
(
elem_type
)
coords
=
self
.
quad_coords
(
elem_type
)
factor
[:]
=
1.
factor
[
coords
[:,
1
]
<
0.5
]
=
.
5
# declares all the parameters that are needed
def
getPushWaveSpeed
(
self
,
params
):
return
np
.
sqrt
((
self
.
lame_lambda
+
2
*
self
.
lame_mu
)
/
self
.
rho
)
# compute small deformation tensor
@staticmethod
def
computeEpsilon
(
grad_u
):
return
0.5
*
(
grad_u
+
np
.
einsum
(
'aij->aji'
,
grad_u
))
# constitutive law
def
computeStress
(
self
,
el_type
,
ghost_type
):
grad_u
=
self
.
getGradU
(
el_type
,
ghost_type
)
sigma
=
self
.
getStress
(
el_type
,
ghost_type
)
n_quads
=
grad_u
.
shape
[
0
]
grad_u
=
grad_u
.
reshape
((
n_quads
,
2
,
2
))
factor
=
self
.
getInternalReal
(
'factor'
)(
el_type
,
ghost_type
)
.
reshape
(
n_quads
)
epsilon
=
self
.
computeEpsilon
(
grad_u
)
sigma
=
sigma
.
reshape
((
n_quads
,
2
,
2
))
trace
=
np
.
einsum
(
'aii->a'
,
grad_u
)
sigma
[:,
:,
:]
=
(
np
.
einsum
(
'a,ij->aij'
,
trace
,
self
.
lame_lambda
*
np
.
eye
(
2
))
+
2.
*
self
.
lame_mu
*
epsilon
)
sigma
[:,
:,
:]
=
np
.
einsum
(
'aij, a->aij'
,
sigma
,
factor
)
# constitutive law tangent modulii
def
computeTangentModuli
(
self
,
el_type
,
tangent_matrix
,
ghost_type
):
n_quads
=
tangent_matrix
.
shape
[
0
]
tangent
=
tangent_matrix
.
reshape
(
n_quads
,
3
,
3
)
factor
=
self
.
getInternalReal
(
'factor'
)(
el_type
,
ghost_type
)
.
reshape
(
n_quads
)
Miiii
=
self
.
lame_lambda
+
2
*
self
.
lame_mu
Miijj
=
self
.
lame_lambda
Mijij
=
self
.
lame_mu
tangent
[:,
0
,
0
]
=
Miiii
tangent
[:,
1
,
1
]
=
Miiii
tangent
[:,
0
,
1
]
=
Miijj
tangent
[:,
1
,
0
]
=
Miijj
tangent
[:,
2
,
2
]
=
Mijij
tangent
[:,
:,
:]
=
np
.
einsum
(
'aij, a->aij'
,
tangent
,
factor
)
# computes the energy density
def
computePotentialEnergy
(
self
,
el_type
):
sigma
=
self
.
getStress
(
el_type
)
grad_u
=
self
.
getGradU
(
el_type
)
nquads
=
sigma
.
shape
[
0
]
stress
=
sigma
.
reshape
(
nquads
,
2
,
2
)
grad_u
=
grad_u
.
reshape
((
nquads
,
2
,
2
))
epsilon
=
self
.
computeEpsilon
(
grad_u
)
energy_density
=
self
.
getPotentialEnergy
(
el_type
)
energy_density
[:,
0
]
=
0.5
*
np
.
einsum
(
'aij,aij->a'
,
stress
,
epsilon
)
# ------------------------------------------------------------------------------
# applies manually the boundary conditions
def
applyBC
(
model
):
nbNodes
=
model
.
getMesh
()
.
getNbNodes
()
position
=
model
.
getMesh
()
.
getNodes
()
displacement
=
model
.
getDisplacement
()
blocked_dofs
=
model
.
getBlockedDOFs
()
width
=
1.
height
=
1.
epsilon
=
1e-8
for
node
in
range
(
0
,
nbNodes
):
if
((
np
.
abs
(
position
[
node
,
0
])
<
epsilon
)
or
(
np
.
abs
(
position
[
node
,
0
]
-
width
)
<
epsilon
)):
blocked_dofs
[
node
,
0
]
=
True
displacement
[
node
,
0
]
=
0
*
position
[
node
,
0
]
+
0.
if
np
.
abs
(
position
[
node
,
1
])
<
epsilon
:
# lower side
blocked_dofs
[
node
,
1
]
=
True
displacement
[
node
,
1
]
=
-
1.
if
np
.
abs
(
position
[
node
,
1
]
-
height
)
<
epsilon
:
# upper side
blocked_dofs
[
node
,
1
]
=
True
displacement
[
node
,
1
]
=
1.
# register the material to the material factory
def
allocator
(
dim
,
option
,
model
,
id
):
return
LocalElastic
(
model
,
id
)
mat_factory
=
aka
.
MaterialFactory
.
getInstance
()
mat_factory
.
registerAllocator
(
"local_elastic"
,
allocator
)
# main parameters
spatial_dimension
=
2
mesh_file
=
'square.msh'
# read mesh
mesh
=
aka
.
Mesh
(
spatial_dimension
)
mesh
.
read
(
mesh_file
)
# parse input file
aka
.
parseInput
(
'material.dat'
)
# init the SolidMechanicsModel
model
=
aka
.
SolidMechanicsModel
(
mesh
)
model
.
initFull
(
_analysis_method
=
aka
.
_static
)
# configure the solver
solver
=
model
.
getNonLinearSolver
()
solver
.
set
(
"max_iterations"
,
2
)
solver
.
set
(
"threshold"
,
1e-3
)
solver
.
set
(
"convergence_type"
,
aka
.
SolveConvergenceCriteria
.
solution
)
# prepare the dumper
model
.
setBaseName
(
"bimaterial"
)
model
.
addDumpFieldVector
(
"displacement"
)
model
.
addDumpFieldVector
(
"internal_force"
)
model
.
addDumpFieldVector
(
"external_force"
)
model
.
addDumpField
(
"strain"
)
model
.
addDumpField
(
"stress"
)
model
.
addDumpField
(
"factor"
)
model
.
addDumpField
(
"blocked_dofs"
)
# Boundary conditions
applyBC
(
model
)
# solve the problem
model
.
solveStep
()
# dump paraview files
model
.
dump
()
epot
=
model
.
getEnergy
(
'potential'
)
print
(
'Potential energy: '
+
str
(
epot
))
Event Timeline
Log In to Comment