Ok, I extended the BCs. I'm getting indices out of bounds error. I suppose this is because I enforce the BCs in DG0, and then solve in CG space. How can I map it from one to another?

Thanks,
Tomasz




Code:

from firedrake import *
import numpy as np

nx = 8; ny = 4; 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) ) * dx
L = Constant(0.) * v * dx

bc = DirichletBC( V, Expression( "sin(pi*x[1]/Ly)", Ly=Ly ), 1 )

class MyBC(DirichletBC):
    def __init__(self, V, value, markers):
        # Call superclass init
        # We provide a dummy subdomain id.
        super(MyBC, self).__init__(V, value, 0)
        # Override the "nodes" property which says where the boundary
        # condition is to be applied.
        self.nodes = np.unique(np.where(markers.dat.data_ro_with_halos == 0)[0])

# Now I can use MyBC to create a "boundary condition" to zero out all
# the nodes that are *not* in the fluid domain:
BC_exclude_beyond_fluid = MyBC( V, 0, I_f )

solve( a==L, phi, bcs=[bc, BC_exclude_beyond_fluid] )



Error listing:

Traceback (most recent call last):
  File "laplace_extended_domain.py", line 51, in <module>
    solve( a==L, phi, bcs=[bc, BC_exclude_beyond_fluid] )
  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 181, in solve
    bc.apply(self._problem.u)
  File "<decorator-gen-297>", line 2, in apply
  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/bcs.py", line 233, in apply
    r.assign(self.function_arg, subset=self.node_set)
  File "/home/tommy/programs/firedrake/local/lib/python2.7/site-packages/pyop2/utils.py", line 64, in __get__
    obj.__dict__[self.__name__] = result = self.fget(obj)
  File "/home/tommy/programs/firedrake/local/lib/python2.7/site-packages/firedrake/bcs.py", line 162, in node_set
    return op2.Subset(self._function_space.node_set, self.nodes)
  File "/home/tommy/programs/firedrake/local/lib/python2.7/site-packages/pyop2/backends.py", line 118, in __call__
    return t(*args, **kwargs)
  File "<decorator-gen-7>", line 2, in __init__
  File "/home/tommy/programs/firedrake/local/lib/python2.7/site-packages/pyop2/utils.py", line 130, in wrapper
    return f(*args, **kwargs)
  File "/home/tommy/programs/firedrake/local/lib/python2.7/site-packages/pyop2/base.py", line 840, in __init__
    (self._indices[0], self._indices[-1], self._superset.total_size))
pyop2.exceptions.SubsetIndexOutOfBounds: Out of bounds indices in Subset construction: [20, 63) not [0, 45)




________________________________________
From: firedrake-bounces@imperial.ac.uk <firedrake-bounces@imperial.ac.uk> on behalf of Lawrence Mitchell <lawrence.mitchell@imperial.ac.uk>
Sent: 15 April 2016 11:22
To: firedrake@imperial.ac.uk
Subject: Re: [firedrake] Laplace equation on subdomain

> On 15 Apr 2016, at 11:08, Tomasz Salwa [RPG] <mmtjs@leeds.ac.uk> wrote:
>
> 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.


The problem is that you've knocked out a part of your operator so that it's:

laplace 0
0       0

Particularly, it's zero on the diagonal outside the subdomain.

To correct this, you'll need to define your own custom boundary condition class that inherits from DirichletBC that identifies the nodes that are outside the subdomain and impose strong conditions there.  I think we discussed this when you visited imperial a while ago.

That way you'll end up with:

laplace 0
0       I

with the identity on the diagonal outside the subdomain instead.

As to your comment (how to write 0*v*dx more succinctly) use:

L = Constant(0)*v*dx

Thanks,

Lawrence

_______________________________________________
firedrake mailing list
firedrake@imperial.ac.uk
https://mailman.ic.ac.uk/mailman/listinfo/firedrake