First of all, your coordinates make no sense, as
will never match with
Please pay special attention to details before posting code.
Your code is also way to complicated, introducing boundary meshes, sub meshes of boundary meshes etc.
Here is a minimal working example of what you want to achieve:
from dolfin import *
import numpy as np
mesh = UnitSquareMesh(4, 4)
VT = FunctionSpace(mesh, "CG", 1)
T_data = np.array([0, 50, 100, 150, 0])
T_coordinates = np.array([[0, 0], [0., 0.25], [0., 0.5], [0., 0.75], [0., 1.]])
assert len(T_data) == len(T_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)
vertices = vt.where_equal(1)
vertex_to_dof_map = vertex_to_dof_map(VT)
T_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(T_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):
T_load.vector()[boundary_dofs[row]] = T_data[col]
print(row, col, boundary_coords[row], T_coordinates[col], T_data[col])
bcT = [DirichletBC(VT, T_load, leftBoundary),
]
u = Function(VT)
[bc.apply(u.vector()) for bc in bcT]
File("mwe.pvd") << u
giving
0 0 [0. 0.] [0. 0.] 0
1 1 [0. 0.25] [0. 0.25] 50
2 2 [0. 0.5] [0. 0.5] 100
3 3 [0. 0.75] [0. 0.75] 150
4 4 [0. 1.] [0. 1.] 0
and function
