# LSQ BC

col_to_index = {}
index = 0
indices = []
total_size = 0

for b in dirichlet_boundaries:
    e = b.element_id()
    
    ref_pts = sample(b)
    total_size+=ref_pts.shape[0]
    
    for l_b in e.bases():
        if is_boundary_node(l_b.global_index):
            if not l_b.global_index in col_to_index:
                col_to_index[l_b.global_index] = index++
                indices.append(index)
                

ii = []
jj = []
vv = []
rhs = zeros((total_size, dim))
g_index = 0

for b in dirichlet_boundaries:
    e = b.element_id()
    ref_pts = sample(b)
    pts = b.eval_gmapping(ref_pts)
    
    bc = eval_dirichlet(pts)
    rhs[g_index:g_index+bc.shape[0], :] = bc
    
    for l_b in e.bases():
        # not boudanry
        if not l_b.global_index in col_to_index:
            continue
        
        b_eval = l_b.eval(ref_pts)
        
        index = col_to_index(l_b.global_index)
        
        for i in range(len(b_eval)):
            ii.appennd(g_index + i)
            index.appennd(g_index)
            vv.append(b_eval[i])
    
    g_index += bc.shape[0]
    
    
M = sparse(ii, jj, vv)
sols = solve(M.T @ M, M.T @ rhs)

for i in range(sols.shape[0]):
    for d in range(dim):
        real_rhs[indices[i]*dim + d] = sols[i, d]