Splitting expressions involving boundary normals

On 25 October 2012 13:21, Øyvind Evju
<email address hidden> wrote:
> New question #212274 on DOLFIN:
> https://answers.launchpad.net/dolfin/+question/212274
>
> I am trying to split an Expression into subcomponents, where this Expression is to be used as a boundary condition.
>
> A simple example is found below. When splitting, I get the error:
> "...
> cell.normal(ufc_cell.local_facet)
> TypeError: expected positive 'int' for argument 2"
>
> Tried to add a try/except block within the eval-method as well, but do not get any results as anticipated. Any ideas?
>
>
> from dolfin import *
>
> class MyExpression(Expression):
> def __init__(self, mesh):
> Expression.__init__(self)
> self.mesh = mesh
>
> def eval_cell(self, value, x, ufc_cell):
> # Find normal
> cell = Cell(self.mesh, ufc_cell.index)
> normal = cell.normal(ufc_cell.local_facet)
>
> value[0] = 0.1*normal.x()
> value[1] = 0.2*normal.y()
> value[2] = 0.3*normal.z()
>
> def value_shape(self):
> return (3,)
>
> class Boundary(SubDomain):
> def inside(self, x, on_boundary):
> return on_boundary
>
> N = 8
> mesh = UnitCube(N,N,N)
> exp = MyExpression(mesh)
>
> # This works fine!
> V = VectorFunctionSpace(mesh, "CG", 1)
> u = Function(V)
> bc = DirichletBC(V, exp, Boundary())
>
> bc.apply(u.vector())
> plot(u, interactive=True)
>
> # Splitting the expression does not work...
> Q = FunctionSpace(mesh, "CG", 1)

Not what you ask for, but this line:
> u = [Function(Q)]*3

is the equivalent of

f = Function(Q)
u = [f, f, f]

not
u = [Function(Q), Function(Q), Function(Q)]
as you probably intended.

You can also do
u = [Function(Q) for _ in range(3)]

About your question, what happens is that exp[i] is an ufl expression
that is projected onto Q, and during that projection exp is called in
the interior of the domain where there is no facet normal available.
You'll have to write a scalar version of MyExpression instead.

Martin

Thank for the reply. I would hate to write a scalar version, so I made a workaround. This only works for flat boundaries though, but that's good enough for me. It is also a lot faster if the Expression is called many times (e.g. in a time loop).

For more general boundaries, I guess you would have to do something more clever.

from dolfin import *

class NormalExpression(Expression):
    def __init__(self, mesh):
        Expression.__init__(self)
        self.mesh = mesh

    def eval_cell(self, value, x, ufc_cell):
        # Find normal
        cell = Cell(self.mesh, ufc_cell.index)
        normal = cell.normal(ufc_cell.local_facet)

        value[0] = normal.x()
        value[1] = normal.y()
        value[2] = normal.z()

    def value_shape(self):
        return (3,)

class MyExpression(Expression):
    def __init__(self, mesh, indicator):
        Expression.__init__(self)
        self.mesh = mesh
        self.ind = indicator
        self.find_average_normal()

    def find_average_normal(self):
        V = VectorFunctionSpace(self.mesh, "CR", 1)
        bc = DirichletBC(V, NormalExpression(self.mesh), self.ind)
        u = Function(V)
        bc.apply(u.vector())

        u0,u1,u2 = u.split(deepcopy=True)

        x = u0.vector().sum()
        y = u1.vector().sum()
        z = u2.vector().sum()

        xyz_len = sqrt(x**2+y**2+z**2)
        if xyz_len == 0:
            self.nx = 0
            self.ny = 0
            self.nz = 0
            return
        self.nx = x/xyz_len
        self.ny = y/xyz_len
        self.nz = z/xyz_len

    def eval_cell(self, value, x, ufc_cell):
        # Find normal
        value[0] = 0.1*self.nx
        value[1] = 0.2*self.ny
        value[2] = 0.3*self.nz

    def value_shape(self):
        return (3,)

class Boundary(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and x[0] < DOLFIN_EPS

N = 8
mesh = UnitCube(N,N,N)
exp = MyExpression(mesh, Boundary())

# This works fine!
V = VectorFunctionSpace(mesh, "CG", 1)
u = Function(V)

bc = DirichletBC(V, exp, Boundary())
bc.apply(u.vector())

plot(u, interactive=True)

# Splitting the expression now works...
Q = FunctionSpace(mesh, "CG", 1)
u = [Function(Q) for _ in range(3)]

for i in range(3):
    bc = DirichletBC(Q, exp[i], Boundary())
    bc.apply(u[i].vector())
    plot(u[i], interactive=True)