Page MenuHomec4science
Files F85316151

darmadi2.py
No OneTemporary
Actions

File Metadata

Created
Sat, Sep 28, 06:46

darmadi2.py
View Options

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Mon Dec 13 15:24:42 2021
@author: ekinkubilay
"""
import numpy as np
import matplotlib.pyplot as plt
from fenics import *
import os
#read the mesh
mesh = Mesh('../Part_geometry/mesh/layer_004.xml')
boundary_mesh = BoundaryMesh(mesh, 'exterior')
velocity = np.array([1, 0, 0])
velocity_mag = np.linalg.norm(velocity)
t = 0.1 #time
alpha = 1
upper = velocity.dot(velocity)*t/(4*alpha) #time
upper_limit = np.min([upper, 5]) #time
lower_limit = 1e-12
V = velocity/(2*alpha)
N = 20
#points = boundary_mesh.coordinates()
#source = np.zeros(np.shape(boundary_mesh.coordinates()))
source_loc = np.array([-1.8, -1.8, 4])
#source[:] = np.array(source_loc)
#disp = source - points
#R_squared = (disp*disp).sum(1)
#R = np.sqrt(R_squared)
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))
return integral_value*(upper-lower)/(3*N)
"""
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_mag*delta.dot(delta))/(4*w*2*alpha))-delta.dot(velocity)/(2*alpha))/(w**1.5)
def flux_function(w, delta):
return np.exp(-w-((velocity_mag*delta.dot(delta))/(4*w*2*alpha))-delta.dot(velocity)/(2*alpha))/(w**2.5)
def get_boundary_normals(mesh):
n = FacetNormal(mesh)
V = VectorFunctionSpace(mesh, "DG", 0)
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
flux_in_cell = {}
flux_in_node = {}
boundary_normals = get_boundary_normals(mesh)
#delta_T = np.zeros(len(points))
"""
for i,j in enumerate(points):
delta_T[i] = integrate(temperature_function, disp[i], lower_limit, upper_limit)
"""
for cell in cells(boundary_mesh):
midpoint = cell.midpoint()[:]
disp = source_loc - midpoint
normal = boundary_normals(midpoint)/np.linalg.norm(boundary_normals(midpoint))
integral = (-1/(4*alpha**2))*velocity.dot(velocity)*integrate(flux_function, disp, lower_limit, upper_limit)
vector = np.multiply(midpoint, normal)*integral
area = cell.volume()
flux_in_cell[cell.index()] = area*vector.dot(normal)*integral
for v in 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()]*(t/3)
#print(v.index(), c.index(), heat_incr[v.index()])
for v in vertices(mesh):
delta_T = integrate(temperature_function, v.point()[:]-source_loc, lower_limit, upper_limit)
#np.savetxt("temp_increase.txt", delta_T)
"""
start = time.time()
flux_in_cell = {}
flux_in_node = {}
boundary_normals = get_boundary_normals(mesh)
for cell in cells(boundary_mesh):
midpoint = cell.midpoint()[:]
disp = np.array([1,1,1]) - midpoint
R_squared = disp.dot(disp)
R = np.sqrt(R_squared)
u = R*V
normal = boundary_normals(midpoint)/np.linalg.norm(boundary_normals(midpoint))
integral = integrate(flux_function, u, R, lower_limit, upper_limit)
vector = np.multiply(midpoint, normal)*integral
area = cell.volume()
#flux_in_cell[cell.index()] = area*vector.dot(normal)*integral
for v in vertices(cell):
try:
flux_in_node[v.index()] += area*vector.dot(normal)*integral
except KeyError:
flux_in_node[v.index()] = flux_in_node.get(v.index(),0)
flux_in_node[v.index()] += area*vector.dot(normal)*integral
print('it took:', time.time()-start)
"""

Event Timeline

c4science · Help