Hi all.

I’m working on a system of non-linear PDEs in a sphere. Accordingly, I have defined everything in spherical coordinates. Since the problem has rotational symmetry, things depend only on r and \theta and the mesh is simply a [0,R]\times[0,\pi] rectangle.

Following the Setting multiple Dirichlet, Neumann, and Robin conditions chapter of the FEniCSx tutorial, I define a custom measure

```
domain = mesh.create_rectangle(MPI.COMM_WORLD, np.array([[0,0],[R,np.pi]]),[10,int(10*np.pi)])
x = ufl.SpatialCoordinate(domain)
boundaries = [(1, lambda x: np.isclose(x[0], 0)), # origin
(2, lambda x: np.isclose(x[0], R)), # surface of the sphere
(3, lambda x: np.isclose(x[1], 0)),
(4, lambda x: np.isclose(x[1], np.pi))]
facet_indices, facet_markers = [], []
fdim = domain.topology.dim - 1
for (marker, locator) in boundaries:
facets = mesh.locate_entities(domain, fdim, locator)
facet_indices.append(facets)
facet_markers.append(np.full_like(facets, marker))
facet_indices = np.hstack(facet_indices).astype(np.int32)
facet_markers = np.hstack(facet_markers).astype(np.int32)
sorted_facets = np.argsort(facet_indices)
facet_tag = mesh.meshtags(domain, fdim, facet_indices[sorted_facets], facet_markers[sorted_facets])
ds = ufl.Measure("ds", domain=domain, subdomain_data=facet_tag)
```

so I can refer only to the r=R part of the boundary when writing the surface contribution to the energy:

```
F = fB(m,x)*Jacob(x)*ufl.dx + fS(m,m0)*Jacob(x)*ufl.ds(2)
```

I then compute the variation of this functional, which gives rise to Robin BCs in that region of the boundary, and solve the non-linear problem:

```
dF = ufl.derivative(F, m, phi) # variation of the energy
JF = ufl.derivative(dF, m, dm) # Jacobian of the energy
problem = fem.petsc.NonlinearProblem(dF, m, J=JF)
```

This procedure, however, doesn’t work:

```
Number of interations: 2
2022-08-09 14:55:33.806 ( 0.110s) [main ] petsc.cpp:677 INFO| PETSc Krylov solver starting to solve system.
2022-08-09 14:55:33.839 ( 0.143s) [main ] petsc.cpp:677 INFO| PETSc Krylov solver starting to solve system.
2022-08-09 14:55:33.846 ( 0.151s) [main ] NewtonSolver.cpp:36 INFO| Newton iteration 2: r (abs) = 0 (tol = 1e-10) r (rel) = -nan(tol = 1e-06)
2022-08-09 14:55:33.847 ( 0.151s) [main ] NewtonSolver.cpp:255 INFO| Newton solver finished in 2 iterations and 2 linear solver iterations.
```

Weirdly enough, if I define the energy by integrating over the complete boundary, `F = fB(m,x)*Jacob(x)*ufl.dx + fS(m,m0)*Jacob(x)*ufl.ds`

, the solver converges well and the result is perfectly compatible with what is expected.

Am I defining the custom measure wrongly? Adding one more layer of complexity, my next step is to define different conditions for r=R,\,\theta<\pi/2 and r=R,\,\theta\geq\pi/2, so I will definitely need to separate the boundary in subregions.