Homec4science
Diffusion LibMultiScale (master) Browse

libmultiscale/test/coupling/LAMMPS_AKANTU/test_qc_static_prestress/mesh/GenerateMesh.py
3fe696d8c6c2master

Last Change
Raw File

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
x_stretch = 1e-2 # one percent of stretch in x direction
alpha = 1+x_stretch
def setParams():
Params = namedtuple("Params", ["a", "latt",
"Lx","Ly","Lz","pad_size"])
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)
Lx = 2
Ly = 8
Lz = 1
pad_size = 3.6
return Params(latt_params, lattice, Lx, Ly, Lz, pad_size)
def make_auxiliary(params, domain):
refinement_point = PointRef3d()
def setter(x,y,z, refinement):
refinement_point[0] = x
refinement_point[1] = y
refinement_point[2] = z
domain.setPointRefinement(refinement_point, refinement)
Lx = params.Lx * params.a.x
Ly = params.Ly * params.a.y
Lp = params.pad_size * params.a.y
Lz = params.Lz * params.a.z
setter(-Lx,Ly,-Lz, 2)
setter( Lx,Ly,-Lz, 2)
setter( Lx,Ly, Lz, 2)
setter(-Lx,Ly, Lz, 2)
setter(-Lx,Lp,-Lz, 0)
setter( Lx,Lp,-Lz, 0)
setter( Lx,Lp, Lz, 0)
setter(-Lx,Lp, Lz, 0)
setter(-Lx,-Ly,-Lz,0)
setter( Lx,-Ly,-Lz,0)
setter( Lx,-Ly, Lz,0)
setter(-Lx,-Ly, Lz,0)
def make_geometries(params):
Lx = params.a.x * params.Lx
Ly = params.a.y * params.Ly
Lz = params.a.z * params.Lz
Lp = params.a.y * params.pad_size
cube = [-Lx, Lx, -Ly, Ly, -Lz, Lz]
overall = Block3d(cube, "overall")
refined_cube = [-Lx, Lx, -Ly, Lp, -Lz, Lz]
refined = Block3d(refined_cube,"refined")
atomic_cube = [-Lx, Lx, -Ly, params.a.y * 0.01, -Lz, Lz]
atomistic = Block3d(atomic_cube,"atomistic")
atomic_skinned_cube = [-Lx, Lx, -Ly, params.a.y * -0.001, -Lz, Lz]
atomistic_skinned = Block3d(atomic_skinned_cube,"atomistic_skinned")
delta = .1
fixed_boundary_cube = [-(Lx+delta), Lx+delta,
Ly-delta, Ly+delta,
-(Lz+delta), Lz+delta]
fixed_boundary = Block3d(fixed_boundary_cube, "fixed_boundary")
return overall, refined, atomistic, atomistic_skinned, fixed_boundary
def make_domain(params, geoms):
overall, refined, atomistic, atomistic_skinned, fixed_boundary= geoms
boundary_special_treatment = True
global center
center = PointRef3d(vecReal([0, 0, 0]))
tolerance = -.001
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()))
print("Domain created with overall geometry: {0}".format(overall.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, geoms):
overall, refined, atomistic, atomistic_skinned, fixed_boundary = geoms
bounds = vecReal(overall.bounds())
for i in [0,2]:
bounds[2*i] += .1
bounds[2*i+1] += .1
#prestresses:
if i == 0:
bounds[2*i] *= alpha
bounds[2*i+1] *= alpha
atoms = PointContainer3d(domain.getAtoms())
print("Writing atoms")
atoms.dump_lammps("atoms.dat",bounds)
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, fixed_boundary = geoms
def print_geom(name, geom, fh):
bound_vals = list(geom.bounds())
bound_vals[0] *= alpha
bound_vals [1] *= alpha
print("GEOMETRY {0} CUBE BBOX {1}".format(
name,
" ".join("{0}".format(val) for val in bound_vals)), file = fh)
geoms = [("non_coupling_atoms", atomistic_skinned),
("all_atoms", overall),
("fixed_boundary", fixed_boundary)]
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(True, False, True)
domain.complement_periodically(-1, 0)
domain.complement_periodically(-1, 2)
domain.getPoints().reportDuplicates(1)
output_mesh(params, domain)
write_geometries(params, geoms)
output_atoms(params, domain, geoms)
if __name__ == "__main__":
print("starting")
main()
print("DONE")
c4science · Help