Thank you， dokken.I will try your method as soon as possible and hope that you can provide modification suggestions afterwards.
Dear dokken，I tried to extract the non-uniform heat flux density boundary conditions to solve the temperature distribution in another area. The code is as follows：
from dolfin import *
import numpy as np
mesh = UnitSquareMesh(10, 10)
VT = FunctionSpace(mesh, "CG", 1)
T_ = TestFunction(VT)
dT = TrialFunction(VT)
q_data = np.linspace(0, 1, 11)
new_y_coordinates = np.linspace(0, 1, 11) #assuming this is the heat flux density extracted from another boundary nodes
q_coordinates = np.column_stack((np.zeros(11), new_y_coordinates))
assert len(q_data) == len(q_coordinates)
class LeftBoundary(SubDomain):
def inside(self, x, on_boundary):
return near(x[0], 0) and on_boundary
leftBoundary = LeftBoundary()
vt = MeshFunction("size_t", mesh, 0, 0)
leftBoundary.mark(vt, 1)
def right(x, on_boundary):
return near(x[0], 1) and on_boundary
vertices = vt.where_equal(1)
vertex_to_dof_map = vertex_to_dof_map(VT)
q_load = Function(VT)
# Get all boundary coordinates
dof_coords = VT.tabulate_dof_coordinates()
boundary_dofs = [vertex_to_dof_map[vertex] for vertex in vertices]
boundary_coords = np.array([dof_coords[boundary_dof] for boundary_dof in boundary_dofs])
# Compute distances between boundary coordinates and T coordinates
points_A = np.expand_dims(boundary_coords, 1)
points_B = np.expand_dims(q_coordinates, 0)
distances = np.sum(np.square(points_A-points_B), axis=2)
is_close = distances < 1e2 * DOLFIN_EPS
positions = np.nonzero(is_close)
for row, col in zip(*positions):
q_load.vector()[boundary_dofs[row]] = q_data[col]
print(row, col, boundary_coords[row], q_coordinates[col], q_data[col])
vt_ = MeshFunction("size_t", mesh, mesh.topology().dim()-1, 0)
leftBoundary.mark(vt_, 1)
ds = Measure('ds', domain=mesh, subdomain_data=vt_)
Delta_T = Function(VT)
aT = dot(grad(dT), grad(T_))*dx
LT = -q_load*T_*ds(1)
bcT = DirichletBC(VT, Constant(1000), right)
solve(aT == LT, Delta_T, bcT)
File("T.pvd") << Delta_T
# Compute heat flux
Q = VectorFunctionSpace(mesh, 'CG', 1)
w = TrialFunction(Q)
v = TestFunction(Q)
aQ = inner(w, v)*dx
LQ = inner(-grad(Delta_T), v)*dx
flux = Function(Q)
solve(aQ == LQ, flux)
File("q.pvd") << flux
The matching of heat flux density with nodes is as follows (looks fine):
0 0 [0. 0.] [0. 0.] 0.0
1 1 [0. 0.1] [0. 0.1] 0.1
2 2 [0. 0.2] [0. 0.2] 0.2
3 3 [0. 0.3] [0. 0.3] 0.30000000000000004
4 4 [0. 0.4] [0. 0.4] 0.4
5 5 [0. 0.5] [0. 0.5] 0.5
6 6 [0. 0.6] [0. 0.6] 0.6000000000000001
7 7 [0. 0.7] [0. 0.7] 0.7000000000000001
8 8 [0. 0.8] [0. 0.8] 0.8
9 9 [0. 0.9] [0. 0.9] 0.9
10 10 [0. 1.] [0. 1.] 1.0
Then I calculated the heat flux density field, and post-processed to obtain the normal heat flux density on the left boundary. Upon comparison, I found that it differs significantly from the set value (especially the red font):
I don’t know if this method is correct. If it’s incorrect, please help me modify it so that the calculated result of heat flux density more accurately matches the set value（How to correctly extract the heat flux density results from a certain boundary and then apply them to solve the temperature distribution in another area). If I am correct, could you explain the reasons for the differences? Thanks for your reply！
1 Like