Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F88866039
GenerateMesh.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
Mon, Oct 21, 01:53
Size
10 KB
Mime Type
text/x-python
Expires
Wed, Oct 23, 01:53 (1 d, 23 h)
Engine
blob
Format
Raw Data
Handle
21831048
Attached To
rLIBMULTISCALE LibMultiScale
GenerateMesh.py
View Options
#!/usr/bin/env python
from
__future__
import
print_function
,
division
import
numpy
as
np
from
collections
import
namedtuple
from
scipy
import
linalg
import
math
from
CADDMesher
import
*
center
=
0
def
setParams
():
Params
=
namedtuple
(
"Params"
,
[
"a"
,
"latt"
,
"height"
,
"width"
,
"thickness"
,
"skinned_limit"
,
"atomistics_limit"
,
"pad_limit"
])
LattParams
=
namedtuple
(
"LattParams"
,
[
"x"
,
"y"
,
"z"
,
"a"
])
a
=
4.04
mil
=
vecReal
([
-
1.
,
1.
,
0.
,
1.
,
1.
,
1.
,
1.
,
1.
,
-
2
])
def
compute_norms
(
matrix
):
def
norm
(
matrix
,
index
):
return
np
.
sqrt
(
matrix
[
index
]
**
2
+
matrix
[
index
+
1
]
**
2
+
matrix
[
index
+
2
]
**
2
)
norms
=
list
()
norms
.
append
(
norm
(
matrix
,
0
))
norms
.
append
(
norm
(
matrix
,
3
))
norms
.
append
(
norm
(
matrix
,
6
))
return
norms
norms
=
compute_norms
(
mil
)
latt_params
=
LattParams
(
a
*
norms
[
0
],
a
*
norms
[
1
],
a
*
norms
[
2
],
a
)
lattice
=
FccLattice3d
(
vecReal
([
a
,
a
,
a
]),
mil
)
height
=
40
*
latt_params
.
y
width
=
18.5
*
latt_params
.
x
thickness
=
6
*
latt_params
.
z
skinned_limit
=
2.9
*
latt_params
.
y
atomistics_limit
=
3
*
latt_params
.
y
pad_thickness
=
2
*
latt_params
.
y
pad_limit
=
atomistics_limit
+
pad_thickness
return
Params
(
latt_params
,
lattice
,
height
,
width
,
thickness
,
skinned_limit
,
atomistics_limit
,
pad_limit
)
def
make_auxiliary
(
params
,
domain
,
upper
):
x_coords
=
-
params
.
a
.
x
,
params
.
width
+
params
.
a
.
x
z_coords
=
-
params
.
a
.
z
,
params
.
thickness
+
params
.
a
.
z
refinement_point
=
PointRef3d
()
sign
=
1
if
upper
else
-
1
for
x
in
x_coords
:
refinement_point
[
0
]
=
x
for
z
in
z_coords
:
refinement_point
[
2
]
=
z
## slip plane:
refinement_point
[
1
]
=
-
params
.
a
.
y
*
sign
refinement
=
0
domain
.
setPointRefinement
(
refinement_point
,
refinement
)
## end of pad:
refinement_point
[
1
]
=
sign
*
(
params
.
pad_limit
-
2
*
params
.
a
.
y
)
refinement
=
0
domain
.
setPointRefinement
(
refinement_point
,
refinement
)
refinement_point
[
1
]
=
sign
*
(
params
.
pad_limit
+
6
*
params
.
a
.
y
)
# this way, the pad limit is also the limit of refined meshes
refinement
=
2.9
domain
.
setPointRefinement
(
refinement_point
,
refinement
)
## some midpoint
refinement_point
[
1
]
=
refinement_point
.
getComponent
(
1
)
+
sign
*
12
*
params
.
a
.
y
domain
.
setPointRefinement
(
refinement_point
,
12
)
## outer edge:
refinement_point
[
1
]
=
sign
*
(
params
.
height
+
params
.
a
.
y
)
domain
.
setPointRefinement
(
refinement_point
,
12
)
def
make_domain
(
params
,
upper
):
if
upper
:
min_max
=
[
-.
5
*
params
.
a
.
x
,
params
.
width
,
0
,
params
.
height
,
0
,
params
.
thickness
]
sign
=
1
y_index
=
3
else
:
min_max
=
[
0.
*
params
.
a
.
x
,
params
.
width
,
-
params
.
height
,
0
,
0
,
params
.
thickness
]
sign
=
-
1
y_index
=
2
## domain geometries
overall
=
Block3d
(
min_max
,
"overall"
)
min_max
[
y_index
]
=
sign
*
params
.
pad_limit
refined
=
Block3d
(
min_max
,
"refined"
)
min_max
[
y_index
]
=
sign
*
params
.
atomistics_limit
atomistic
=
Block3d
(
min_max
,
"atomistic"
)
min_max
[
y_index
]
=
sign
*
params
.
skinned_limit
atomistic_skinned
=
Block3d
(
min_max
,
"atomistic_skinned"
)
boundary_special_treatment
=
True
global
center
center
=
PointRef3d
(
vecReal
([
params
.
a
.
x
*.
05
,
params
.
a
.
y
*.
01
,
params
.
a
.
z
*.
05
]))
center
=
PointRef3d
(
vecReal
([
0
,
0
,
0
]))
tolerance
=
.
1
domain
=
Domain3d
(
overall
,
refined
,
atomistic
,
atomistic_skinned
,
params
.
latt
,
boundary_special_treatment
,
center
,
.
1
)
make_auxiliary
(
params
,
domain
,
upper
)
return
domain
def
output_mesh
(
params
,
domain
,
upper
):
mesh
=
domain
.
getMesh
(
100
)
print
(
"Quality for {0} mesh: q_min = {1}, q_max = {2}"
.
format
(
"upper"
if
upper
else
"lower"
,
mesh
.
getRadiusRatioMin
(),
mesh
.
getRadiusRatioMax
()))
ratiovec
=
mesh
.
computeRadiusRatios
()
print
(
"ratiovecsize = {0}"
.
format
(
len
(
ratiovec
)))
ratios
=
np
.
array
(
ratiovec
)
print
(
"ratios: max = {0}, min = {1}"
.
format
(
ratios
.
max
(),
ratios
.
min
()))
print
(
"Writing paraview"
)
mesh
.
dump_paraview
(
"{0}_domain"
.
format
(
"upper"
if
upper
else
"lower"
))
print
(
"Writing gmsh"
)
boundary_direction
=
1
boundary_val
=
params
.
height
if
upper
else
-
params
.
height
mesh
.
tagBoundaryPlane
(
boundary_direction
,
boundary_val
)
mesh
.
dump_msh
(
"{0}_domain_akantu.msh"
.
format
(
"upper"
if
upper
else
"lower"
))
def
getLatticeCoords
(
params
,
atom
,
center
):
A
=
.
5
*
np
.
matrix
([
np
.
array
(
params
.
latt
.
getLatticeStdVector
(
0
))
*
params
.
latt
.
getConstant
(
0
),
np
.
array
(
params
.
latt
.
getLatticeStdVector
(
1
))
*
params
.
latt
.
getConstant
(
1
),
np
.
array
(
params
.
latt
.
getLatticeStdVector
(
2
))
*
params
.
latt
.
getConstant
(
2
)])
.
T
b
=
np
.
array
([
atom
[
i
]
-
center
[
i
]
for
i
in
range
(
3
)])
latt_coordinates
=
linalg
.
solve
(
A
,
b
)
rounded_coordinates
=
np
.
array
([
round
(
coord
)
for
coord
in
latt_coordinates
])
error
=
np
.
abs
(
latt_coordinates
-
rounded_coordinates
)
.
sum
()
if
error
>
1e-8
:
print
(
"Error = {0}"
.
format
(
error
))
return
latt_coordinates
def
output_atoms
(
params
,
domain_upper
,
domain_lower
):
def
getFilteredAtoms
(
domain
):
atoms
=
PointContainer3d
(
domain
.
getAtoms
())
bounds
=
atoms
.
computeBounds
()
atom_geom
=
Block3d
(
bounds
,
"atomic geometry"
)
# fix for periodicity: in fem the left (min) node is the master of its
# periodic equivalent at the right (max)
## temp ##atom_geom.shift(PointRef3d(vecReal([params.a.x/4, params.a.y/4, 0*-params.a.z*.1])))
## temp ##atoms.filter(atom_geom)
## temp ##bounds = list(bounds)
print
(
bounds
)
return
atoms
,
atom_geom
upper_atoms
,
upper_geom
=
getFilteredAtoms
(
domain_upper
)
upper_atoms
.
dump_lammps
(
"upper_atoms.dat"
)
upper_atoms
.
dump_paraview
(
"upper_atoms"
)
##min_y = 10
##for i in range(upper_atoms.getSize()):
## if upper_atoms.getPoint(i)[1] < min_y:
## min_y = upper_atoms.getPoint(i)[1]
## getLatticeCoords(params, upper_atoms.getPoint(i), center)
##
lower_atoms
,
lower_geom
=
getFilteredAtoms
(
domain_lower
)
lower_atoms
.
dump_lammps
(
"lower_atoms.dat"
)
lower_atoms
.
dump_paraview
(
"lower_atoms"
)
max_y
=
-
10
##for i in range(lower_atoms.getSize()):
## if lower_atoms.getPoint(i)[1] > max_y:
## max_y = lower_atoms.getPoint(i)[1]
## getLatticeCoords(params, lower_atoms.getPoint(i), upper_atoms.getPoint(12))
##
##print("min_y = {0}, max_y = {1}".format(min_y, max_y))
atoms
=
PointContainer3d
(
upper_atoms
)
## now add lower atoms to atoms
for
point_id
in
range
(
lower_atoms
.
getSize
()):
atoms
.
addPoint
(
lower_atoms
.
getPoint
(
point_id
))
## update the bounds
bounds
=
vecReal
(
6
)
for
i
in
range
(
3
):
bounds
[
2
*
i
+
0
]
=
min
(
upper_geom
.
bounds
()[
2
*
i
+
0
],
lower_geom
.
bounds
()[
2
*
i
+
0
])
bounds
[
2
*
i
+
1
]
=
max
(
upper_geom
.
bounds
()[
2
*
i
+
1
],
lower_geom
.
bounds
()[
2
*
i
+
1
])
## expand them so they're the next higher multiple of a lattice (or lattice + .5 burgers)
addition
=
[
-.
248
,
0
,
0
]
for
i
in
range
(
3
):
size
=
bounds
[
2
*
i
+
1
]
-
bounds
[
2
*
i
]
rounded_size
=
(
math
.
ceil
(
size
/
params
.
a
[
i
]
-
addition
[
i
])
+
addition
[
i
])
*
params
.
a
[
i
]
delta
=
rounded_size
-
size
bounds
[
2
*
i
+
0
]
-=
.
5
*
delta
bounds
[
2
*
i
+
1
]
+=
.
5
*
delta
print
((
"old size in {0} direction = {2:.3f} ({3:.3f} lattices, with lattice "
"size in {0} = {4:.3f})
\n
"
"new size in {0} direction = {5:.3f} ({6} lattices, with lattice "
"size in {0} = {4:.3f})
\n
"
)
.
format
((
'x'
,
'y'
,
'z'
)[
i
],
i
,
size
,
size
/
params
.
a
[
i
],
params
.
a
[
i
],
rounded_size
,
rounded_size
/
params
.
a
[
i
]))
print
(
list
(
bounds
))
print
(
"size of upper container = {0}, lower container = {1}, total container = {2}"
.
format
(
upper_atoms
.
getSize
(),
lower_atoms
.
getSize
(),
atoms
.
getSize
()))
print
(
"Writing atoms"
)
atoms
.
dump_paraview
(
"all_atoms"
)
atoms
.
dump_lammps
(
"all_atoms.dat"
,
bounds
)
##FIX THE GEOMETRY
print
(
"Writing interface and pad atoms"
)
print
(
params
)
upper_interface_atoms
=
PointContainer3d
(
"upper_interface_atoms"
)
upper_pad_atoms
=
PointContainer3d
(
"upper_pad_atoms"
)
## The automatic detection of interface atoms detects also atoms as on the
## lateral interface as interface atoms, even though the are pad atoms due
## to the periodic boundary conditions. Therefore, we explicitely create an
## interface bounding box:
min_max
=
[
-.
5
*
params
.
a
.
x
,
params
.
width
,
params
.
atomistics_limit
-.
5
,
params
.
atomistics_limit
+.
5
,
0
,
params
.
thickness
]
upper_interface_box
=
Block3d
(
min_max
)
domain_upper
.
compute_coupling_atoms
(
upper_interface_atoms
,
upper_pad_atoms
,
upper_interface_box
)
upper_interface_atoms
.
dump_paraview
(
"upper_interface_atoms"
)
upper_pad_atoms
.
dump_paraview
(
"upper_pad_atoms"
)
upper_interface_atoms
.
dump_text
(
"upper_interface_atoms.txt"
)
upper_pad_atoms
.
dump_text
(
"upper_pad_atoms.txt"
)
lower_interface_atoms
=
PointContainer3d
(
"lower_interface_atoms"
)
lower_pad_atoms
=
PointContainer3d
(
"lower_pad_atoms"
)
## The automatic detection of interface atoms detects also atoms as on the
## lateral interface as interface atoms, even though the are pad atoms due
## to the preriodic boundary conditions. Therefore, we explicitely create an
## interface bounding box:
min_max
=
[
0.
*
params
.
a
.
x
,
params
.
width
,
-
params
.
atomistics_limit
-.
5
,
-
params
.
atomistics_limit
+.
5
,
0
,
params
.
thickness
]
lower_interface_box
=
Block3d
(
min_max
)
domain_lower
.
compute_coupling_atoms
(
lower_interface_atoms
,
lower_pad_atoms
,
lower_interface_box
)
lower_interface_atoms
.
dump_paraview
(
"lower_interface_atoms"
)
lower_pad_atoms
.
dump_paraview
(
"lower_pad_atoms"
)
lower_interface_atoms
.
dump_text
(
"lower_interface_atoms.txt"
)
lower_pad_atoms
.
dump_text
(
"lower_pad_atoms.txt"
)
def
main
():
params
=
setParams
()
def
getDomain
(
upper
):
domain
=
make_domain
(
params
,
upper
)
domain
.
fill_points
(
True
,
False
,
True
)
domain
.
complement_periodically
(
1
,
0
)
domain
.
complement_periodically
(
1
,
2
)
domain
.
getPoints
()
.
reportDuplicates
(
1
)
output_mesh
(
params
,
domain
,
upper
)
return
domain
output_atoms
(
params
,
getDomain
(
upper
=
True
),
getDomain
(
upper
=
False
))
print
(
getLatticeCoords
(
params
,
PointRef3d
(
vecReal
([
4.5707382335898199
,
-
32.584956372739398
,
59.045766746369502
])),
center
))
if
__name__
==
"__main__"
:
print
(
"starting"
)
main
()
print
(
"DONE"
)
Event Timeline
Log In to Comment