Homec4science
Diffusion LibMultiScale (features/cadd) Browse

libmultiscale/test/coupling/MD1D_MECA1D/test_quasi_continuum_stability/mesh/cadd1d.py
c8dd2d0ba119features/cadd

Last Change
Raw File

cadd1d.py
View Options

#!/usr/bin/env python
import numpy as np
import numpy.matlib as mp
import numpy.linalg as alg
from scipy.optimize import minimize, fmin_bfgs
#import matplotlib.pyplot as plt
import sys
class CADD1D(object):
def __init__(self, stiffness, mass, lattice, nb_fem, nb_md, nb_pad):
"""
represents a 1D cadd chain with the fem at the right end and attached at the fem
"""
# self.nb_interaction = 2 # hardcoded 2nd neighbour interaction for md
self.nb_fem = nb_fem + nb_pad
self.nb_md = nb_md
self.nb_pad = nb_pad
self.nb_points = nb_fem + nb_md + nb_pad
self.stiffness = float(stiffness)
self.mass = float(mass)
self.lattice = float(lattice)
# Fem is formulated in displacements, md in positions
self.__points = np.arange(0, (self.nb_points)*self.lattice, self.lattice)
self.position = self.points.copy()
self.generateConnectivity()
self.generateLennardJones()
@property
def points(self):
return self.__points
@property
def displacement(self):
return self.position - self.points
def generateLennardJones(self):
nb_interaction = 2 # hardcoded 2nd neighbour interaction for md
def zeta(x):
return sum([1./i**x for i in range(1, nb_interaction+1)])
self.sigma = self.lattice*(2*zeta(12)/zeta(6))**(-1./6)
self.epsilon = self.stiffness*self.lattice/(24*(-7*self.sigma**6*zeta(6)/self.lattice**7
+ 26* self.sigma**12*zeta(12)/self.lattice**13))
def repr_material_params(self, linestart = ""):
outstr = list()
outstr.append("input parameters:")
outstr.append(" stiffness k = {0}".format(self.stiffness))
outstr.append(" mass m = {0}".format(self.mass))
outstr.append(" lattice param = {0}".format(self.lattice))
outstr.append(" LJ sigma = {0}".format(self.sigma))
outstr.append(" LJ epsilon = {0}".format(self.epsilon))
outstr.append("derived parameters:")
outstr.append(" stable timestep = {0}".format(np.sqrt(self.mass/self.stiffness)))
return "{0}{1}".format(linestart, "\n{0}".format(linestart).join(outstr))
def MDforce(self, p1, p2):
r = (self.position[p1] - self.position[p2])#not taking the abs takes care of the sign of the force
r6 = r**6
r7 = r6*r
r13 = r6*r7
return (24*self.epsilon*(self.sig6zeta6/r7 - 2*self.sig12zeta12/r13))
def MDenergy(self, p1, p2):
r6 = (self.position[p1] - self.position[p2])**6
r12 = r6**2
return 4*self.epsilon*(self.sig12/r12 - self.sig6/r6)
def FEMforce(self, p1, p2):
delta_u = (self.position[p1] - self.position[p2]) - (self.points[p1] - self.points[p2])
return -delta_u*self.stiffness
def FEMenergy(self, p1, p2):
delta_u = (self.position[p1] - self.position[p2]) - (self.points[p1] - self.points[p2])
return delta_u**2*self.stiffness/2.
def generateConnectivity(self, ):
""" The connectivity here is the dual of what's usually meant: a neighbourlist
"""
self.neighs = list()
def clean(iterable):
return tuple([index for index in iterable if 0 <= index < self.nb_points])
for pt_id in range(self.nb_fem):
self.neighs.append(clean([pt_id - 1, pt_id + 1]))
def isFemMaster(self, index):
return 0 <= index < self.nb_fem
def isFem(self, index):
return index <= self.nb_fem
def isMDMaster(self, index):
return not (self.isFemMaster(index) )
def isMD(self, index):
return index >= self.nb_fem - self.nb_pad
def isPad(self, index):
return self.isMD(index) and self.isFemMaster(index)
def isInterface(self, index):
return self.isMDMaster(index) and self.isFem(index)
def isCoupled(self, index):
return self.isMD(index) and self.isFem(index)
def dump(self, filename_pattern):
def sub_dump(name, filter_fun):
name = "{0}-{1}.txt".format(filename_pattern, name)
with open(name, "w") as fh:
print("# file: {0}".format(name), file=fh)
print(self.repr_material_params("# "), file=fh)
for pt in [self.points[index] for index in range(self.nb_points) if filter_fun(index)]:
print(pt, file=fh)
sub_dump("interface_atoms", self.isInterface)
sub_dump("pad_atoms", self.isPad)
sub_dump("all_atoms", self.isMD)
sub_dump("nodes", self.isFem)
# compute the geoms
upper_lim = max(self.points) + self.lattice/2.
lower_lim = min(self.points) - self.lattice/2.
refin_lim = min([self.points[id] for id in range(self.nb_points)
if self.isPad(id)]) - self.lattice/2.
atoms_lim = self.points[self.nb_fem] - self.lattice/2.
skinn_lim = atoms_lim + self.lattice
geoms = list()
geoms.append(("overall", (lower_lim, upper_lim)))
geoms.append(("refined", (refin_lim, upper_lim)))
geoms.append(("atomistic", (atoms_lim, upper_lim)))
geoms.append(("skinned", (skinn_lim, upper_lim)))
with open("{0}-geometries.config".format(filename_pattern), "w") as fh:
print("Section GEOMETRIES AtomsUnits", file=fh)
print("LET stableTimestep = {0}".format(np.sqrt(self.mass/self.stiffness)), file=fh)
print("GEOMETRY full CUBE BBOX -10000 10000", file=fh)
for geom in geoms:
print("GEOMETRY {0} CUBE BBOX {1[0]} {1[1]}".format(*geom), file=fh)
print("endSection", file=fh)
def main():
cadd_args = dict()
cadd_args["stiffness"] = 10
cadd_args["mass"] = 2
cadd_args["lattice"] = 4
cadd_args["nb_fem"] = 10
cadd_args["nb_md"] = 10
cadd_args["nb_pad"] = 10
cadd = CADD1D(**cadd_args)
cadd.dump("test")
if __name__ == "__main__":
main()
c4science · Help