Page MenuHomec4science
Files F96153367

flux_calculation.py
No OneTemporary
Actions

File Metadata

Created
Mon, Dec 23, 05:36

flux_calculation.py
View Options

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Tue Dec 14 13:44:10 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')
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
# for cell in cells(boundary_mesh):
# print('cell: ', cell)
# print('cell index: ', cell.index())
# print('cell midpoint:', cell.midpoint()[:])
# print('vertex coordinates: ', cell.get_vertex_coordinates())
# print('volume:', cell.volume())
# print('normal:', nh(cell.midpoint()[:]))
# heat_incr = {}
# dummy_vector = np.array([1,1,1])
# for v in vertices(boundary_mesh):
# for c in cells(v):
# normal = nh(c.midpoint()[:])/np.linalg.norm(nh(c.midpoint()[:]))
# area = c.volume()
# incr = (1/3)*area*dummy_vector.dot(normal)
# heat_incr[v.index()] = heat_incr.get(v.index(),0) + incr
# #print(v.index(), c.index(), heat_incr[v.index()])
dummy_vector = np.array([1,1,1])
flux_in_cell = {}
flux_in_node = {}
for cell in cells(boundary_mesh):
midpoint = cell.midpoint()[:]
normal = nh(midpoint)/np.linalg.norm(nh(midpoint))
#normal = nh(cell.midpoint()[:])/np.linalg.norm(nh(cell.midpoint()[:]))
area = cell.volume()
incr = area*dummy_vector.dot(normal)
flux_in_cell[cell.index()] = (1/3)*area*dummy_vector.dot(normal)
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()]/3
#print(v.index(), c.index(), heat_incr[v.index()])

Event Timeline

c4science · Help