Page MenuHomec4science
Files F88846184

GenerateMesh.py
No OneTemporary
Actions

File Metadata

Created
Sun, Oct 20, 23:33

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

c4science · Help