Dear Firedrakers,

As an intermediate part of fluid-structure interaction problem, I would like to solve Laplace equation (with zero Neumann boundary conditions) on a subdomain of the full domain. For that reason I define the step function that is 1 in the subdomain of interest and 0 elsewhere. Multiplying LHS of the weak form (with trial function) by the step function causes a solver error. How can I correct it? Code and error listing are below.

Thank you,

Tomasz

from firedrake import *

nx = 17; ny = 8; Lx = 2.; Ly = 1.

mesh = RectangleMesh( nx, ny, Lx, Ly )
V = FunctionSpace( mesh, "CG", 1 )
V_DG0 = FunctionSpace( mesh, "DG", 0 )

def Heaviside( zero_on_left=True, Lp=0. ):
    sign = 1.
    if zero_on_left == False:
        sign = -1.
    Hx = Function( V_DG0 )
    Hx_expr = Expression( "0.5*(1.0 + copysign( 1.0, sign*(x[0]-Lp) ))",
                            sign=sign, Lp=Lp )
    Hx.interpolate(Hx_expr)
    return Hx;

I_f = Heaviside( zero_on_left=False, Lp=1. )

phi = Function( V )
trial = TrialFunction( V )
v = TestFunction( V )

a = dot( grad(trial), grad(v) ) * I_f * dx   # step function causes error here !!!
L = Function(V) * v * dx   # BTW can I put zero here in a simpler manner?

bc = DirichletBC( V, Expression( "sin(pi*x[1]/Ly)", Ly=Ly ), 1 )
solve( a==L, phi, bcs=bc )

Error listing:

...

Traceback (most recent call last):
  File "laplace_extended_domain.py", line 36, in <module>
    solve( a==L, phi, bcs=bc )
  File "/home/tommy/programs/firedrake/local/lib/python2.7/site-packages/firedrake/solving.py", line 120, in solve
    _solve_varproblem(*args, **kwargs)
  File "/home/tommy/programs/firedrake/local/lib/python2.7/site-packages/firedrake/solving.py", line 147, in _solve_varproblem
    solver.solve()
  File "<decorator-gen-295>", line 2, in solve
  File "/home/tommy/programs/firedrake/local/lib/python2.7/site-packages/pyop2/profiling.py", line 203, in wrapper
    return f(*args, **kwargs)
  File "/home/tommy/programs/firedrake/local/lib/python2.7/site-packages/firedrake/variational_solver.py", line 190, in solve
    solving_utils.check_snes_convergence(self.snes)
  File "/home/tommy/programs/firedrake/local/lib/python2.7/site-packages/firedrake/solving_utils.py", line 62, in check_snes_convergence
    %s""" % (snes.getIterationNumber(), msg))
RuntimeError: Nonlinear solve failed to converge after 0 nonlinear iterations.
Reason:
   Inner linear solve failed to converge after 0 iterations with reason: unknown reason (petsc4py enum incomplete?)