Defining Multiple Neumann BCs

[KQEFEM]

Hi all,
I am struggling with defining multiple Neumann BCs for my geometry.

I want to choose the inner and outer boundary such that they have different conditions. Here is the model I am considering

$$ 
\begin{cases}
- \nabla \cdot ( \kappa \nabla u) = f \quad &\text{on } \Omega \\

\kappa(\nabla u \cdot \underline{n} )= -q( u - T_1) \quad &\text{on } \partial \Omega_1 \\

\kappa(\nabla u \cdot \underline{n} )= -q( u - T_2) \quad &\text{on } \partial \Omega_2

\end{cases}, 
$$ 

where q is a given constant.

This has the following following bilinear form:
$$ \sum_\mathbf{K}\kappa \int_{\mathbf{K}} \nabla u \cdot \nabla v , dx - \kappa \sum_{\mathbf{K}, i }\int_{\partial \Omega_{\mathbf{K}i}} v (-q(u-T_i) , ds = \sum\mathbf{K}\int_{\mathbf{K}} f v , dx.
$$

I tried loading in my mesh and defining the boundary basis separately but this gives me some funky results!

  # Precompute global basis at quadrature points
  basis = fem.Basis(mesh, element_type)  # shorthand for CellBasis
  # Identify the inner and outer boundaries
  inner_boundary = (
      meshconfig["boundary_indices_dict"]["inner_boundary"] - 1
  )  # mesh.boundaries["inner"]
  outer_boundary = (
      meshconfig["boundary_indices_dict"]["outer_boundary"] - 1
  )  # mesh.boundaries["outer"]

  # Create separate boundary bases for inner and outer boundaries
  boundary_basis_inner = basis.boundary(inner_boundary)
  boundary_basis_outer = basis.boundary(outer_boundary)

# functions etc

  # Assemble the stiffness matrix (A) and load vector (b)
  A = laplace.assemble(basis)
  b = rhs.assemble(basis)

  neumann_mat_inner = boundary_integral_neumann_inner.assemble(boundary_basis_inner)
  neumann_mat_outer = boundary_integral_neumann_outer.assemble(boundary_basis_outer)

  neumann_rhs_inner = boundary_integral_neumann_rhs_inner.assemble(
      boundary_basis_inner
  )
  neumann_rhs_outer = boundary_integral_neumann_rhs_outer.assemble(
      boundary_basis_outer
  )

  b += neumann_rhs_inner + neumann_rhs_outer
  A += neumann_mat_inner + neumann_mat_outer

A = LHS +/- neumman_mat
b = RHS +/- neumman_rhs

X = fem.solve(A,b) 

In your other examples (in the discussion) you use condense but I do not see how this takes into account the BC as my understanding was this reduces the problem to a Dirichlet?

Some advise on how to properly implement this would be ideal!