Hi, I am trying to solve the following linear equation in a simple 2D domain (like a unit square):

\nabla ( u \nabla T) = 0 \ in \ \Omega, \\ T + u (n, \nabla T) = 0 \ on \ \Gamma_1 \\ \mathrm{(no \ bcs \ on \ } \Gamma_2 == \partial \Omega - \Gamma_1 )

Here `u`

is the unknown scalar field, and `T`

is a known scalar field (not sure that is relevant, but I obtained `T`

as a solution to a different problem).

So far I could only write something like this for the weak formulation:

\int_\Omega - u (\nabla v, \nabla T) dx + \int_{\Gamma_2} v u (n, \nabla T) ds = \int_{\Gamma_2} v T ds

which translates to the code:

```
from dolfinx.mesh import create_unit_square
from mpi4py import MPI
import numpy as np
from dolfinx import mesh as dolf_mesh
from dolfinx import fem
import ufl
from ufl import dx, grad, inner
from dolfinx.fem.petsc import LinearProblem
mesh = create_unit_square(MPI.COMM_WORLD, 100, 100)
V_e = fem.functionspace(mesh, ("Lagrange", 2))
# An example field `T`
T = fem.Function(V_e)
T.interpolate(lambda x: ((x[0] - 0.5)**2.0 + (x[1] - 1.5) ** 2.0)**0.5)
def bot_boundary(x):
return np.isclose(x[1], 0.0)
def rest_boundary(x):
return np.logical_not(bot_boundary(x))
# Got this from https://github.com/FEniCS/dolfinx/blob/1ccffb0aa545634ea6e160236a7b043dc2a7fdd2/python/test/unit/fem/test_assemble_domains.py#L121
bottom_facets = dolf_mesh.locate_entities_boundary(
mesh, mesh.topology.dim - 1, bot_boundary
)
bottom_vals = np.full(bottom_facets.shape, 1, np.intc)
rest_facets = dolf_mesh.locate_entities_boundary(
mesh, mesh.topology.dim - 1, rest_boundary
)
rest_vals = np.full(rest_facets.shape, 2, np.intc)
indices = np.hstack((bottom_facets, rest_facets))
values = np.hstack((bottom_vals, rest_vals))
indices, pos = np.unique(indices, return_index=True)
marker = dolf_mesh.meshtags(mesh, mesh.topology.dim - 1, indices, values[pos])
ds = ufl.Measure("ds", subdomain_data=marker, domain=mesh)
v = ufl.TestFunction(V_e)
u = ufl.TrialFunction(V_e)
a = -u * inner(grad(v), grad(T)) * dx
a += v * u * inner(grad(T), ufl.FacetNormal(mesh)) * ds(2)
L = v * T * ds(1)
problem = LinearProblem(a, L, bcs=[])
uh = problem.solve()
```

However the results just look rubbish, and in fact the solution sometimes consists of `inf`

s depending on my mesh (with `petsc_options={"ksp_type": "preonly", "pc_type": "lu"}`

in the problem definition).

I am really struggling to understand what I am doing wrong and would appreciate any tips of how to solve this. Thank you!