Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F88846184
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
Sun, Oct 20, 23:33
Size
5 KB
Mime Type
text/x-python
Expires
Tue, Oct 22, 23:33 (1 d, 22 h)
Engine
blob
Format
Raw Data
Handle
21826411
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"
,
"cubesize"
,
"skinned_limit"
,
"atomistics_limit"
,
"pad_limit"
])
LattParams
=
namedtuple
(
"LattParams"
,
[
"x"
,
"y"
,
"z"
,
"a"
])
a
=
4.0452598
# minimised value computed in lammps for the mendelev potential
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
)
cubesize
=
8
atomistics_limit
=
2.1
skinned_limit
=
1.9
pad_thickness
=
1.8
pad_limit
=
atomistics_limit
+
pad_thickness
return
Params
(
latt_params
,
lattice
,
cubesize
,
skinned_limit
,
atomistics_limit
,
pad_limit
)
def
make_auxiliary
(
params
,
domain
):
coords
=
np
.
array
([
-
1.
,
1.
],
dtype
=
float
)
refinement_point
=
PointRef3d
()
for
x_val
in
coords
:
for
y_val
in
coords
:
for
z_val
in
coords
:
x
=
x_val
*
params
.
a
.
x
y
=
y_val
*
params
.
a
.
y
z
=
z_val
*
params
.
a
.
z
def
setter
(
factor
,
refinement
):
refinement_point
[
0
]
=
factor
*
x
refinement_point
[
1
]
=
factor
*
y
refinement_point
[
2
]
=
factor
*
z
domain
.
setPointRefinement
(
refinement_point
,
refinement
)
setter
(
1.01
*
params
.
cubesize
,
4
)
setter
(
params
.
atomistics_limit
,
0
)
def
make_geometries
(
params
):
min_max
=
[
-
params
.
a
.
x
,
params
.
a
.
x
,
-
params
.
a
.
y
,
params
.
a
.
y
,
-
params
.
a
.
z
,
params
.
a
.
z
]
overall
=
Block3d
([
params
.
cubesize
*
val
for
val
in
min_max
],
"overall"
)
refined
=
Block3d
([
params
.
pad_limit
*
val
for
val
in
min_max
],
"refined"
)
atomistic
=
Block3d
([
params
.
atomistics_limit
*
val
for
val
in
min_max
],
"atomistic"
)
atomistic_skinned
=
Block3d
([
params
.
skinned_limit
*
val
for
val
in
min_max
],
"atomistic_skinned"
)
return
overall
,
refined
,
atomistic
,
atomistic_skinned
def
make_domain
(
params
,
geoms
):
overall
,
refined
,
atomistic
,
atomistic_skinned
=
geoms
boundary_special_treatment
=
True
global
center
center
=
PointRef3d
(
vecReal
([
0
,
0
,
0
]))
tolerance
=
.
1
domain
=
Domain3d
(
overall
,
refined
,
atomistic
,
atomistic_skinned
,
params
.
latt
,
boundary_special_treatment
,
center
,
tolerance
)
print
(
"Domain created with atomistic geometry: {0}"
.
format
(
atomistic_skinned
.
bounds
()))
make_auxiliary
(
params
,
domain
)
return
domain
def
output_mesh
(
params
,
domain
):
mesh
=
domain
.
getMesh
(
100
)
print
(
"Quality for mesh: q_min = {0}, q_max = {1}"
.
format
(
mesh
.
getRadiusRatioMin
(),
mesh
.
getRadiusRatioMax
()))
print
(
"Writing paraview"
)
mesh
.
dump_paraview
(
"domain"
)
print
(
"Writing gmsh"
)
mesh
.
dump_msh
(
"domain_akantu.msh"
)
def
output_atoms
(
params
,
domain
):
atoms
=
PointContainer3d
(
domain
.
getAtoms
())
print
(
"Writing atoms"
)
atoms
.
dump_lammps
(
"atoms.dat"
)
atoms
.
dump_paraview
(
"atoms"
)
print
(
"Writing interface and pad atoms"
)
interface_atoms
=
PointContainer3d
(
"interface_atoms"
)
pad_atoms
=
PointContainer3d
(
"pad_atoms"
)
domain
.
compute_coupling_atoms
(
interface_atoms
,
pad_atoms
)
interface_atoms
.
dump_paraview
(
"interface_atoms"
)
pad_atoms
.
dump_paraview
(
"pad_atoms"
)
interface_atoms
.
dump_text
(
"interface_atoms.txt"
)
pad_atoms
.
dump_text
(
"pad_atoms.txt"
)
def
write_geometries
(
params
,
geoms
):
overall
,
refined
,
atomistic
,
atomistic_skinned
=
geoms
def
print_geom
(
name
,
geom
,
fh
):
print
(
"GEOMETRY {0} CUBE BBOX {1}"
.
format
(
name
,
" "
.
join
(
"{0}"
.
format
(
val
)
for
val
in
geom
.
bounds
())),
file
=
fh
)
geoms
=
[(
"non_coupling_atoms"
,
atomistic_skinned
),
(
"all_atoms"
,
overall
)]
with
open
(
"geometries.config"
,
"w"
)
as
fh
:
print
(
"Section GEOMETRIES MetalUnits"
,
file
=
fh
)
print
(
"GEOMETRY full CUBE BBOX -1000 1000 -1000 1000 -1000 1000"
,
file
=
fh
)
[
print_geom
(
name
,
geom
,
fh
)
for
(
name
,
geom
)
in
geoms
]
print
(
"GEOMETRY coupling_region SUB IDS {0} {1}"
.
format
(
geoms
[
1
][
0
],
geoms
[
0
][
0
]),
file
=
fh
)
print
(
"endSection"
,
file
=
fh
)
def
main
():
params
=
setParams
()
geoms
=
make_geometries
(
params
)
domain
=
make_domain
(
params
,
geoms
)
domain
.
fill_points
(
False
,
False
,
False
)
domain
.
getPoints
()
.
reportDuplicates
(
1
)
# should never report duplicates
output_mesh
(
params
,
domain
)
write_geometries
(
params
,
geoms
)
output_atoms
(
params
,
domain
)
if
__name__
==
"__main__"
:
print
(
"starting"
)
main
()
print
(
"DONE"
)
Event Timeline
Log In to Comment