Page MenuHomec4science
Files F85220636

flux_calculation.py
No OneTemporary
Actions

File Metadata

Created
Fri, Sep 27, 14:46

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')
boundary_mesh.init_cell_orientations(Expression(('x[0]*0.9', 'x[1]*0.9', 'x[2]'), degree=1))
cell_orientations = boundary_mesh.cell_orientations()
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()])
heat_incr1 = {}
heat_incr2 = {}
for cell in cells(boundary_mesh):
normal = nh(cell.midpoint()[:])/np.linalg.norm(nh(cell.midpoint()[:]))
area = cell.volume()
incr = area*dummy_vector.dot(normal)
heat_incr1[cell.index()] = (1/3)*area*dummy_vector.dot(normal)
for v in vertices(boundary_mesh):
for c in cells(v):
heat_incr2[v.index()] = heat_incr2.get(v.index(),0) + heat_incr1[c.index()]/3
#print(v.index(), c.index(), heat_incr[v.index()])

Event Timeline

c4science · Help