Homec4science
Diffusion Additive Manufacturing Work (temperature_mesh_coarsening) Browse

R11910/Thermal_model/helpers.py
3e9a1d30b01ftemperature_mesh_coarse…

Last Change
Raw File

helpers.py
View Options

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Wed Dec 22 16:13:33 2021
@author: kubilay
"""
from fenics import *
import numpy as np
import numexpr as ne
import matplotlib.pyplot as plt
import pyslm
from pyslm import hatching as hatching
def moving_integral_1(delta_square, delta_velocity, velocity_mag, lower, upper, alpha, limits):
N = 20
A_min, B_min, A_max, B_max = limits
delta = np.sqrt(delta_square)
dof = len(delta)
time_ranges = np.zeros((dof, N))
time_ranges = np.linspace(np.amax([A_min*delta+B_min, lower*np.ones(dof)], axis=0), np.amin([A_max*delta+B_max, upper*np.ones(dof)], axis=0), N).T
integrant = np.zeros((N, len(delta_square)))
for i,t in enumerate(time_ranges.T):
integrant[i] = ne.evaluate('exp(-delta_square/(4*alpha*t) - delta_velocity/(2*alpha) - velocity_mag*t/(4*alpha))/(t**1.5)')
integrant *= np.ptp(time_ranges, axis=1)/(3*N)
B = np.ones(N)
B[1:-1:2] = 4
B[2:-1:2] = 2
return B.dot(integrant)
def moving_integral_2(delta_square, delta_velocity, velocity_mag, lower, upper, alpha, limits):
N = 20
A_min, B_min, A_max, B_max = limits
delta = np.sqrt(delta_square)
dof = len(delta)
time_ranges = np.zeros((dof, N))
time_ranges = np.linspace(np.amax([A_min*delta+B_min, lower*np.ones(dof)], axis=0), np.amin([A_max*delta+B_max, upper*np.ones(dof)], axis=0), N).T
integrant = np.zeros((N, len(delta_square)))
for i,t in enumerate(time_ranges.T):
integrant[i] = ne.evaluate('exp(-delta_square/(4*alpha*t) - delta_velocity/(2*alpha) - velocity_mag*t/(4*alpha))/(t**2.5)')
integrant *= np.ptp(time_ranges, axis=1)/(3*N)
B = np.ones(N)
B[1:-1:2] = 4
B[2:-1:2] = 2
return B.dot(integrant)
def get_boundary_normals(mesh):
n = FacetNormal(mesh)
V = VectorFunctionSpace(mesh, "CG", 1)
u = TrialFunction(V)
v = TestFunction(V)
a = inner(u,v)*ds
l = inner(n, v)*ds
A = assemble(a, keep_diagonal=True)
L = assemble(l)
A.ident_zeros()
nh = Function(V)
solve(A, nh.vector(), L)
#File("nh.pvd") << nh
return nh
def analytical_flux(delta_square, delta_velocity, velocity, normal_vectors, delta, lower, t, alpha, limits):
term_integral = moving_integral_1(delta_square, delta_velocity, velocity.dot(velocity), lower, t, alpha, limits)
term_coeff = ne.evaluate('sum(normal_vectors*velocity, axis=1)')
flux_vector = ne.evaluate('term_integral*term_coeff')
term_integral = moving_integral_2(delta_square, delta_velocity, velocity.dot(velocity), lower, t, alpha, limits)
term_coeff = ne.evaluate('sum(normal_vectors*delta, axis=1)')
flux_vector += ne.evaluate('term_integral*term_coeff')
return flux_vector
def run_analytical_fields(mesh_coordinates, layer, t, alpha, coeff):
temperature_vector = np.zeros(len(mesh_coordinates))
for scan in layer:
if t > scan.t_start:
delta_square = np.einsum('ij,ij->i', mesh_coordinates-scan.initial_loc-scan.velocity*(t-scan.t_start), mesh_coordinates-scan.initial_loc-scan.velocity*(t-scan.t_start))
delta_velocity = (mesh_coordinates-scan.initial_loc-scan.velocity*(t-scan.t_start)).dot(scan.velocity)
temperature_vector += moving_integral_1(delta_square, delta_velocity, scan.velocity.dot(scan.velocity), 0.000000001, t-scan.t_start, alpha, scan.limits)
if t > scan.t_end:
temperature_vector -= moving_integral_1(delta_square, delta_velocity, scan.velocity.dot(scan.velocity), 0.000000001, t-scan.t_end, alpha, scan.limits)
return temperature_vector*coeff
def run_correction_fields(mesh_coordinates, normal_vectors, layer, t, alpha, coeff):
flux_vector = np.zeros(len(mesh_coordinates))
for scan in layer:
if t > scan.t_start:
delta = mesh_coordinates-scan.initial_loc-scan.velocity*(t-scan.t_start)
delta_square = np.einsum('ij,ij->i', mesh_coordinates-scan.initial_loc-scan.velocity*(t-scan.t_start), mesh_coordinates-scan.initial_loc-scan.velocity*(t-scan.t_start))
delta_velocity = delta.dot(scan.velocity)
flux_vector += analytical_flux(delta_square, delta_velocity, scan.velocity, normal_vectors, delta, 0.000000001, t-scan.t_start, alpha, scan.limits)
if t > scan.t_end:
flux_vector -= analytical_flux(delta_square, delta_velocity, scan.velocity, normal_vectors, delta, 0.000000001, t-scan.t_end, alpha, scan.limits)
return flux_vector*coeff
def get_interpolation_limits(u_maximum, v, alpha):
def function(w):
return np.exp(-u**2/(4*alpha*w) - w*v**2/(4*alpha))/w**1.5
u_range = np.linspace(0.000001, u_maximum)
maximums = []
minimums = []
u_maxs = []
u_mins = []
for u in u_range:
w = np.linspace(0.000001, 0.05, 10000)
values = function(w)
limit = 0.001*function(w[np.argmax(values)])*np.ones(len(w))
ids = np.argwhere(np.diff(np.diff(np.sign(values-limit))))
"""
plt.plot(w, values, marker='.')
plt.scatter(w[np.argmax(values)], function(w[np.argmax(values)]), color='r', s=50)
plt.scatter(w[ids], values[ids], color='g', s=50)
plt.show()
"""
if len(ids)==1:
maximums.append(w[ids[0]])
u_maxs.append(u)
else:
minimums.append(w[ids[0]])
maximums.append(w[ids[-1]])
u_maxs.append(u)
u_mins.append(u)
A = np.array([[np.sum(np.power(u_maxs,2)), np.sum(u_maxs)],[np.sum(u_maxs), len(u_maxs)]])
B = np.array([np.asarray(u_maxs).dot(np.asarray(maximums))[0],np.sum(maximums) ])
A_max,B_max = np.linalg.inv(A).dot(B)
A = np.array([[np.sum(np.power(u_mins,2)), np.sum(u_mins)],[np.sum(u_mins), len(u_mins)]])
B = np.array([np.asarray(u_mins).dot(np.asarray(minimums))[0],np.sum(minimums) ])
A_min,B_min = np.linalg.inv(A).dot(B)
plt.plot(u_maxs, maximums, label="numerical upper limits")
plt.plot(u_mins, minimums, label="numerical lower limits")
plt.plot(u_range, 1.05*(A_max*u_range + B_max), label="interpolated upper limits")
plt.plot(u_range, 0.95*(A_min*u_range + B_min), label="numerical lower limits")
plt.xlabel("u (m)")
plt.ylabel("t (s)")
plt.legend()
plt.savefig("Figures/moving_integration_limits.png")
plt.close()
return A_min, B_min, A_max, B_max
def get_scan_strategy(layer_no, contour_velocity_mag, infill_velocity_mag, contour_scan_limits, infill_scan_limits, contour_infill_delay):
# Imports the part and sets the geometry to an STL file (frameGuide.stl)
solidPart = pyslm.Part('inversePyramid')
solidPart.setGeometry('/home/ekinkubilay/Documents/AM/Part_geometry/layer/layer_001.stl')
#solidPart.origin[0] = 0
#solidPart.origin[1] = 0
#solidPart.scaleFactor = 1.0
#solidPart.rotation = [0, 0, 0]
#solidPart.dropToPlatform()
# Create a StripeHatcher object for performing any hatching operations
myHatcher = hatching.Hatcher()
# Set the base hatching parameters which are generated within Hatcher
if layer_no%2 == 1:
myHatcher.hatchAngle = 90
elif layer_no%2 == 0:
myHatcher.hatchAngle = 0
myHatcher.spotCompensation = 0.05
myHatcher.numInnerContours = 0
myHatcher.numOuterContours = 1
myHatcher.scanContourFirst = True
myHatcher.hatchDistance = 0.05
myHatcher.hatchSortMethod = hatching.AlternateSort()
# Set the layer thickness
#layerThickness = 0.04 # [mm]
#infill_velocity_mag = 0.8
#contour_velocity_mag = 0.9
#alpha = 8e-6
#contour_infill_delay = get_interpolation_limits(0.00377/2, infill_velocity_mag, alpha)
#contour_scan_limits = get_interpolation_limits(0.00377/2, contour_velocity_mag, alpha)
#Perform the hatching operations
print('Hatching Started')
layers = []
contour_coordinates = []
infill_coordinates = []
geomSlice = solidPart.getVectorSlice(solidPart.boundingBox[5]-1e-7)
# Hatch the boundary using myHatcher
layer = myHatcher.hatch(geomSlice)
# Plot the layer geometries using matplotlib
#pyslm.visualise.plot(layer, plot3D=False, plotArrows=True, plotOrderLine=False, plotColorbar=False)
for i in layer.getContourGeometry():
for j in i.coords:
contour_coordinates.append(np.array([j[0], j[1] , solidPart.boundingBox[5]]))
contour_coordinates = np.array(contour_coordinates)/1000
for i in layer.getHatchGeometry():
for j in i.coords:
infill_coordinates.append(np.array([j[0], j[1] , solidPart.boundingBox[5]]))
infill_coordinates = np.array(infill_coordinates)/1000
t = 0.00001
contour_layer = Layer()
for i in range(len(contour_coordinates)-1):
del_x = contour_coordinates[i+1] - contour_coordinates[i]
vel = contour_velocity_mag*del_x/np.linalg.norm(del_x)
del_t = np.linalg.norm(del_x)/np.linalg.norm(vel)
s = Scan(t, t+del_t, contour_coordinates[i], vel, contour_scan_limits)
contour_layer.addScan(s)
t += del_t
t += contour_infill_delay
infill_layer = Layer()
i = 0
while i <= (len(infill_coordinates)-2):
del_x = infill_coordinates[i+1] - infill_coordinates[i]
vel = infill_velocity_mag*del_x/np.linalg.norm(del_x)
del_t = np.linalg.norm(del_x)/np.linalg.norm(vel)
s = Scan(t, t+del_t, infill_coordinates[i], vel, infill_scan_limits)
infill_layer.addScan(s)
t += del_t
i += 1
try:
del_x = infill_coordinates[i+1] - infill_coordinates[i]
vel = infill_velocity_mag*del_x/np.linalg.norm(del_x)
del_t = np.linalg.norm(del_x)/np.linalg.norm(vel)
t += del_t
i += 1
except IndexError:
break
print('Completed Hatching')
contour_layer.combineLayers(infill_layer)
return contour_layer
class Scan:
def __init__(self, t_start, t_end, initial_loc, velocity, limits):
self.t_start = t_start
self.t_end = t_end
self.initial_loc = initial_loc
self.velocity = velocity
self.limits = limits
class Layer:
def __init__(self):
self._scans = list()
def addScan(self, scan):
self._scans.append(scan)
def removeScan(self, scan):
self._scans.remove(scan)
def numberOfScans(self):
return len(self._scans)
def combineLayers(self, layer_adding):
for s in layer_adding:
self.addScan(s)
def totalLayerTime(self):
t_max = 0
for s in self._scans:
if s.t_end > t_max:
t_max = s.t_end
return t_max
def __iter__(self):
return LayerIterator(self)
class LayerIterator():
def __init__(self, layer):
self._layer = layer
self._index = 0
def __next__(self):
if self._index < len(self._layer._scans):
result = self._layer._scans[self._index]
self._index += 1
return result
raise StopIteration
"""
def integral_1(delta_square, delta_velocity, velocity_mag, lower, upper, alpha):
N = 100
time_range = np.linspace(lower, upper, N)
integrant = np.zeros((N, len(delta_square)))
for i,t in enumerate(time_range):
integrant[i] = ne.evaluate('exp(-delta_square/(4*alpha*t) - delta_velocity/(2*alpha) - velocity_mag*t/(4*alpha))/(t**1.5)')
B = np.ones(N)
B[1:-1:2] = 4
B[2:-1:2] = 2
B *= (upper-lower)/(3*N)
return B.dot(integrant)
def integral_2(delta_square, delta_velocity, velocity_mag, lower, upper, alpha):
N = 100
time_range = np.linspace(lower, upper, N)
integrant = np.zeros((N, len(delta_square)))
for i,t in enumerate(time_range):
integrant[i] = ne.evaluate('exp(-delta_square/(4*alpha*t) - delta_velocity/(2*alpha) - velocity_mag*t/(4*alpha))/(t**2.5)')
B = np.ones(N)
B[1:-1:2] = 4
B[2:-1:2] = 2
B *= (upper-lower)/(3*N)
return B.dot(integrant)
def integral_1_selective_adaptive(delta_square, delta_velocity, velocity_mag, lower, upper, alpha, err):
result = integral_1(delta_square, delta_velocity, velocity_mag, lower, upper, alpha)
quad_1 = integral_1(delta_square, delta_velocity, velocity_mag, lower, upper*0.5, alpha)
quad_2 = integral_1(delta_square, delta_velocity, velocity_mag, upper*0.5, upper, alpha)
#plt.plot(quad_1+quad_2)
indices = np.invert(np.isclose(quad_1+quad_2, result, rtol=err, atol=1e-8))
#print(len(result[indices]))
while not np.all( np.isclose(quad_1[indices]+quad_2[indices], result[indices], rtol=err, atol=1e-8)):
result[indices] = quad_1[indices] + quad_2[indices]
quad_1[indices] = integral_1_selective_adaptive(delta_square[indices], delta_velocity[indices], velocity_mag, lower, upper*0.5, alpha, err)
quad_2[indices] = integral_1_selective_adaptive(delta_square[indices], delta_velocity[indices], velocity_mag, upper*0.5, upper, alpha, err)
#print(result)
return quad_1 + quad_2
def integral_2_selective_adaptive(delta_square, delta_velocity, velocity_mag, lower, upper, alpha, err):
result = integral_2(delta_square, delta_velocity, velocity_mag, lower, upper, alpha)
quad_1 = integral_2(delta_square, delta_velocity, velocity_mag, lower, upper*0.5, alpha)
quad_2 = integral_2(delta_square, delta_velocity, velocity_mag, upper*0.5, upper, alpha)
#plt.plot(quad_1+quad_2)
indices = np.invert(np.isclose(quad_1+quad_2, result, rtol=err, atol=1e-8))
#print(len(result[indices]))
while not np.all( np.isclose(quad_1[indices]+quad_2[indices], result[indices], rtol=err, atol=1e-8)):
result[indices] = quad_1[indices] + quad_2[indices]
quad_1[indices] = integral_2_selective_adaptive(delta_square[indices], delta_velocity[indices], velocity_mag, lower, upper*0.5, alpha, err)
quad_2[indices] = integral_2_selective_adaptive(delta_square[indices], delta_velocity[indices], velocity_mag, upper*0.5, upper, alpha, err)
#print(result)
return quad_1 + quad_2
def integrate(func, delta, lower, upper):
B = np.ones(N)
B[1:-1:2] = 4
B[2:-1:2] = 2
w = np.linspace(lower, upper, N)
integral_value = B.dot(func(w, delta))*(upper-lower)/(3*N)
#print(integral_value*(upper-lower)/(3*N))
return integral_value
def temperature_function(w, u, R):
return np.exp(-w-((u**2)/(4*w))-R*V)/(w**1.5)
def temperature_function(w, delta):
return np.exp(-w-((velocity.dot(velocity)*delta.dot(delta))/(4*w*4*alpha**2))-delta.dot(velocity)/(2*alpha))/(w**1.5)
def temperature_function2(w, delta, nu):
print(w)
print(-delta.dot(delta)/w)
print(-nu.dot(nu)*w/4)
print(-delta.dot(nu)/2)
return np.exp(-delta.dot(delta)/w - nu.dot(nu)*w/4 - delta.dot(nu)/2)/(w**1.5)
def temperature_function3(t, disp, velocity, alpha):
return np.exp(-((disp.dot(disp/(4*alpha)))/t)-disp.dot(velocity)/(2*alpha)-velocity.dot(velocity)*t/(4*alpha))/(t**1.5)
def flux_function(w, delta):
return np.exp(-w-((velocity.dot(velocity)*delta.dot(delta))/(4*w*4*alpha**2))-delta.dot(velocity)/(2*alpha))/(w**2.5)
class MyExpression(UserExpression):
def __init__(self, time, velocity, source_loc, alpha):
super().__init__(time, velocity, source_loc, degree=1)
self.time = time
self.alpha = alpha
self.velocity = velocity
self.upper_limit = np.min([self.velocity.dot(self.velocity)*self.time/(4*self.alpha),5]) #time
self.lower_limit = 0.000001
self.source_loc = source_loc
self.B = np.ones(N)
self.B[1:-1:2] = 4
self.B[2:-1:2] = 2
self.w = np.linspace(self.lower_limit, self.time+self.lower_limit, N)
def eval(self, value, x):
self.disp = (np.array([x[0], x[1], x[2]])-self.source_loc)*1e-3
#print("disp", self.disp)
#self.u = self.disp.dot(self.disp)*self.velocity.dot(self.velocity)/(4*self.alpha**2)
#print(self.u)
#value[0] = 1.4e4*integrate(temperature_function, self.disp, self.lower_limit, self.upper_limit)
value[0] = (self.time-self.lower_limit)/(3*N)*self.B.dot(temperature_function3(self.w, self.disp, self.velocity, self.alpha))
#value[0] = 1.4e4*((self.upper_limit-self.lower_limit)/(3*N))*self.B.dot(np.exp(-self.w-((self.velocity.dot(self.velocity)*self.disp.dot(self.disp))/(4*self.w*4*self.alpha**2))-self.disp.dot(self.velocity)/(2*self.alpha))/(self.w**1.5))
# self.u < 120:
#print(self.u)
# value[0] = 1.4e0*integrate(temperature_function, self.disp, self.lower_limit, self.upper_limit)
#else:
# value[0] = 0
def value_shape(self):
return ()
def analytical_flux(boundary_mesh, boundary_normals, source_loc, velocity, t, alpha, t_offset):
flux_in_cell = {}
for cell in cells(boundary_mesh):
midpoint = cell.midpoint()[:]
normal = boundary_normals(midpoint)/np.linalg.norm(boundary_normals(midpoint))
area = cell.volume()
disp = midpoint - source_loc-velocity*t
flux_in_cell[cell.index()] = normal.dot(integral_2(midpoint-source_loc-velocity*t, 0.000001, t-t_offset, velocity, alpha)*disp
+ integral_1(midpoint-source_loc-velocity*t, 0.000001, t-t_offset, velocity, alpha)*velocity)
flux = np.zeros(len(boundary_mesh.coordinates()[:]))
for i,v in enumerate(vertices(boundary_mesh)):
for c in cells(v):
#flux_in_node[v.index()] = flux_in_node.get(v.index(),0) + flux_in_cell[c.index()]
flux[i] += (1/3)*flux_in_cell[c.index()]
return flux
def integral_1(delta, lower, upper, velocity, alpha):
#if delta.dot(velocity) > 10*alpha:
# return 0
#if upper > 20*alpha/(velocity.dot(velocity)):
# upper = 20*alpha/(velocity.dot(velocity))
N = 40
B = np.ones(N)
B[1:-1:2] = 4
B[2:-1:2] = 2
t = np.linspace(lower, upper, N)
value = np.exp(-((delta.dot(delta/(4*alpha)))/t)-delta.dot(velocity)/(2*alpha)-velocity.dot(velocity)*t/(4*alpha))/(t**1.5)
integral_value = B.dot(value)*(upper-lower)/(3*N)
return integral_value
def integral_2(delta, lower, upper, velocity, alpha):
N = 40
B = np.ones(N)
B[1:-1:2] = 4
B[2:-1:2] = 2
t = np.linspace(lower, upper, N)
value = np.exp(-((delta.dot(delta/(4*alpha)))/t)-delta.dot(velocity)/(2*alpha)-velocity.dot(velocity)*t/(4*alpha))/(t**2.5)
integral_value = B.dot(value)*(upper-lower)/(3*N)
return integral_value
def run_analytical_field(mesh_coordinates, scan, t, alpha, coeff):
delta_square = np.einsum('ij,ij->i', mesh_coordinates-scan.initial_loc-scan.velocity*(t-scan.t_start), mesh_coordinates-scan.initial_loc-scan.velocity*(t-scan.t_start))
delta_velocity = (mesh_coordinates-scan.initial_loc-scan.velocity*(t-scan.t_start)).dot(scan.velocity)
temperature_vector = integral_1(delta_square, delta_velocity, scan.velocity.dot(scan.velocity), scan.t_start, t-scan.t_start, alpha)
if t > scan.t_end:
temperature_vector -= integral_1(delta_square, delta_velocity, scan.velocity.dot(scan.velocity), scan.t_start, t-scan.t_start-scan.t_end, alpha)
return temperature_vector*coeff
def run_correction_field(mesh_coordinates, normal_vectors, scan, t, alpha, coeff):
delta = mesh_coordinates-scan.initial_loc-scan.velocity*(t-scan.t_start)
delta_square = np.einsum('ij,ij->i', mesh_coordinates-scan.initial_loc-scan.velocity*(t-scan.t_start), mesh_coordinates-scan.initial_loc-scan.velocity*(t-scan.t_start))
delta_velocity = delta.dot(scan.velocity)
flux_vector = analytical_flux(delta_square, delta_velocity, scan.velocity, normal_vectors, delta, scan.t_start, t-scan.t_start, alpha)
if t > scan.t_end:
flux_vector -= analytical_flux(delta_square, delta_velocity, scan.velocity, normal_vectors, delta, scan.t_start, t-scan.t_start-scan.t_end, alpha)
return flux_vector*coeff
def integral_1111(delta_square, delta_velocity, velocity_mag, lower, upper, alpha):
N = 80
time_range = np.logspace(np.log10(lower), np.log10(upper), N)
integrant = np.zeros((N, len(delta_square)))
for i,t in enumerate(time_range):
integrant[i] = ne.evaluate('exp(-delta_square/(4*alpha*t) - delta_velocity/(2*alpha) - velocity_mag*t/(4*alpha))/(t**1.5)')
result = np.trapz(integrant, time_range, axis=0)
B = np.ones(N)
B[1:-1:2] = 4
B[2:-1:2] = 2
B *= (upper-lower)/(3*N)
return result
def integral_1(delta_square, delta_velocity, velocity_mag, lower, upper, alpha):
N = 40
time_range = np.linspace(1e-5, 5e-4, N)
integrant = np.zeros((N, len(delta_square)))
for i,t in enumerate(time_range):
integrant[i] = ne.evaluate('exp(-delta_square/(4*alpha*t) - delta_velocity/(2*alpha) - velocity_mag*t/(4*alpha))/(t**1.5)')
time_range = np.linspace(5e-4, upper, N)
integrant = np.zeros((N, len(delta_square)))
for i,t in enumerate(time_range):
integrant[i] += ne.evaluate('exp(-delta_square/(4*alpha*t) - delta_velocity/(2*alpha) - velocity_mag*t/(4*alpha))/(t**1.5)')
B = np.ones(N)
B[1:-1:2] = 4
B[2:-1:2] = 2
B *= (upper-lower)/(3*N)
return B.dot(integrant)
"""
def create_scan_strategy(mesh):
#determine the limits of the domain in m
x_max = np.max(mesh.coordinates()[:,0])
x_min = np.min(mesh.coordinates()[:,0])
y_max = np.max(mesh.coordinates()[:,1])
y_min = np.min(mesh.coordinates()[:,1])
z_max = np.max(mesh.coordinates()[:,2])
z_min = np.min(mesh.coordinates()[:,2])
infill_offset = 1e-4
contour_offset = 1e-4
y_partition = np.arange(y_min+contour_offset, y_max-contour_offset+infill_offset, step=infill_offset)
start_points = np.empty((len(y_partition), 3))
start_points[:,1] = y_partition
start_points[:,0] = x_min + contour_offset
start_points[:,2] = z_max
end_points = np.empty((len(y_partition), 3))
end_points[:,1] = y_partition
end_points[:,0] = x_max - contour_offset
end_points[:,2] = z_max
return zip(start_points, end_points)
c4science · Help