Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F88867145
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, 02:01
Size
13 KB
Mime Type
text/x-python
Expires
Wed, Oct 23, 02:01 (1 d, 23 h)
Engine
blob
Format
Raw Data
Handle
21821701
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
import
sys
import
argparse
import
os
from
CADDMesher
import
*
center
=
0
refinement_factor
=
1
# if we want full refinement, I toggle this to zero
is_fully_md
=
False
def
getArgs
():
parser
=
argparse
.
ArgumentParser
(
description
=
"Generate a CADD mesh for a simple coupled dislocation"
)
parser
.
add_argument
(
"--md_size"
,
dest
=
"md_size"
,
default
=
3
,
type
=
float
,
help
=
"Half size of the MD zone in number of lattices "
"from the slipplane. This will usually be an integer"
)
parser
.
add_argument
(
"--height"
,
dest
=
"height"
,
default
=
40
,
type
=
float
,
help
=
"Half size of the overall domain in number of "
"lattices from the slipplane. This will will usually "
"be an integer"
)
parser
.
add_argument
(
"--thickness"
,
dest
=
"thickness"
,
default
=
6
,
type
=
int
,
help
=
"Thickness of the overall domain in number of "
"lattices. Has to be an integer"
)
parser
.
add_argument
(
"--width"
,
dest
=
"width"
,
default
=
18.5
,
type
=
float
,
help
=
"Width of overall domain in lattices. This has to "
"be an multiple of 1/2"
)
args
=
parser
.
parse_args
()
if
not
args
.
width
*
2
==
float
(
int
(
args
.
width
*
2
)):
parser
.
error
(
"Invalid width: needs to be a multiple of 1/2"
)
return
args
def
setParams
():
args
=
getArgs
()
Params
=
namedtuple
(
"Params"
,
[
"a"
,
"latt"
,
"height"
,
"width"
,
"thickness"
,
"skinned_limit"
,
"atomistics_limit"
,
"pad_limit"
])
LattParams
=
namedtuple
(
"LattParams"
,
[
"x"
,
"y"
,
"z"
,
"a"
])
a
=
4.04527
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
)
int_height
=
args
.
height
int_pad_thickness
=
2
height
=
int_height
*
latt_params
.
y
width
=
args
.
width
*
latt_params
.
x
thickness
=
args
.
thickness
*
latt_params
.
z
pad_thickness
=
int_pad_thickness
*
latt_params
.
y
md_size
=
args
.
md_size
# md_size = 0 means full refinement
if
md_size
==
0
:
is_fully_md
=
True
global
refinement_factor
refinement_factor
=
0
md_size
=
int_height
-
int_pad_thickness
atomistics_limit
=
md_size
*
latt_params
.
y
skinned_limit
=
(
md_size
-.
1
)
*
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
*
refinement_factor
domain
.
setPointRefinement
(
refinement_point
,
refinement
)
## end of pad:
refinement_point
[
1
]
=
sign
*
(
params
.
pad_limit
-
2
*
params
.
a
.
y
)
refinement
=
0
*
refinement_factor
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
*
refinement_factor
domain
.
setPointRefinement
(
refinement_point
,
refinement
)
## some midpoint
refinement_point
[
1
]
=
refinement_point
.
getComponent
(
1
)
+
sign
*
12
*
params
.
a
.
y
refinement
=
12
*
refinement_factor
domain
.
setPointRefinement
(
refinement_point
,
refinement
)
## outer edge:
refinement_point
[
1
]
=
sign
*
(
params
.
height
+
params
.
a
.
y
)
refinement
=
12
*
refinement_factor
domain
.
setPointRefinement
(
refinement_point
,
refinement
)
def
make_domain
(
params
,
upper
):
if
upper
:
min_max
=
[
0.
,
params
.
width
,
0
,
params
.
height
,
0
,
params
.
thickness
]
sign
=
1
y_index
=
3
else
:
min_max
=
[
0.
,
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"
)
if
is_fully_md
:
boundary_special_treatment
=
False
else
:
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
,
tolerance
)
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
if
upper
:
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
=
[
0.
,
0.
,
0.
]
for
i
in
range
(
3
):
size
=
bounds
[
2
*
i
+
1
]
-
bounds
[
2
*
i
]
# The box size can only take a multiple of burgers vectors in x dimension.
# One burgers is 1/2 of a lattice
resolution
=
2.
rounded_size
=
(
math
.
ceil
(
resolution
*
(
size
/
params
.
a
[
i
]
-
addition
[
i
]))
/
resolution
+
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
():
cwd
=
os
.
getcwd
()
input_path
=
"input_files"
if
not
os
.
path
.
isdir
(
input_path
):
os
.
mkdir
(
input_path
)
os
.
chdir
(
input_path
)
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
))
os
.
chdir
(
cwd
)
if
__name__
==
"__main__"
:
print
(
"starting"
)
main
()
print
(
"DONE"
)
Event Timeline
Log In to Comment