Homec4science
Diffusion LibMultiScale (lammps_2020) Browse

libmultiscale/test/coupling/LAMMPS_AKANTU/test_quasi_continuum_parallel/GenerateMesh.py
301d6de84779lammps_2020

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
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"
)
c4science · Help