Function in boundary condition, defined on boundary only

variational formulation
1

Greetings. I have some function. It is included in the right part of the system of nonlinear PDEs under the normal derivative. The problem is that I have a function defined on boundary only. And multiplication on the test function gives me TypeError: '<' not supported between instances of 'Mesh' and 'Mesh'

import itertools

from utilities import Normal

# noinspection PyUnresolvedReferences
from dolfin import dx, ds
from dolfin import *


class Normal(UserExpression):
    def __init__(self, mesh, dimension, **kwargs):
        self.mesh = mesh
        self.dimension = dimension
        super().__init__(**kwargs)

    def eval_cell(self, values, x, ufc_cell):
        values[0] = 0
        values[1] = 0
        for i in range(self.dimension):

            if x[i] < DOLFIN_EPS:
                values[i] = -1
            if x[i] + DOLFIN_EPS > 1:
                values[i] = 1
        for i, j in itertools.combinations(list(range(self.dimension)), 2):
            if abs(values[i]) == abs(values[j]) == 1:
                values[0] *= 0.5
                values[1] *= 0.5

    def value_shape(self):
        return (2,) if self.dimension == 2 else (3,)


omega = UnitSquareMesh(10, 10)
omega_b = BoundaryMesh(omega, 'exterior')
finite_element = FiniteElement("Lagrange", omega.ufl_cell(), 2)

state_space = FunctionSpace(omega, finite_element * finite_element)
simple_space = FunctionSpace(omega, finite_element)
vector_space = VectorFunctionSpace(omega, 'Lagrange', 2)
boundary_vector_space = VectorFunctionSpace(omega_b, 'Lagrange', 1)
boundary_simple_space = FunctionSpace(omega_b, 'Lagrange', 1)
normal = interpolate(Normal(omega, dimension=2), boundary_vector_space)
normal_2 = FacetNormal(omega_b)
theta = project(
    Expression('pow((x[0]-0.5), 2) - 0.5*x[1] + 0.75', element=simple_space.ufl_element()),
    simple_space
)
grad_t_b = interpolate(project(grad(theta), vector_space), boundary_vector_space)

theta_n = project(inner(grad_t_b, normal), boundary_simple_space)
v, h = TestFunctions(state_space)
state = Function(state_space)
a, b = split(state)

a_eq = inner(grad(a), grad(v)) * dx + inner(a ** 4 - b, v) * dx
b_eq = inner(grad(b), grad(h)) * dx + inner(b - a ** 4, h) * dx
a_src = 0.3 * v * ds
b_src = theta_n * h * ds
solve(a_eq + b_eq - a_src - b_src == 0, state)

What am I missing? How can I include boundary-only function into the boundary conditions?

Thanks in advance.

2

You cannot mix meshes in a single assembly, one of the reasons is that the integration measure is not well defined in that case.
Could you please explain why you cannot use

n = FacetNormal(omega)
theta_n = inner(grad(theta),n)
b_src = theta_n * h * ds
1 Like
3

Thanks for your reply! I will try your suggestion.
To your question - I have a few reasons:

  1. Issue not with the normal derivative itself but rather with “function defined on boundary only” - on my actual case it’s a 3-rd type boundary condition. So I have to use a combination of the normal derivative plus some other function. Yeah, I guess I could refactor my code a bit and maybe get a result. I will do and report what happened. :slight_smile:

  2. I would like to be able to manipulate freely with my functions to trace what’s going on. FacetNormal gives Integral of type cell cannot contain a ReferenceNormal outside of the problem formulation scope (I’m sorry if I’m wrong here) so I personally want to avoid to use it.

4

I have previously posted a recipe on how to project the FacetNormal to an appropriate volume function space, see: How to plot normal unit vector of faces in a 2D mesh? - #2 by dokken
This function would only give you a normal on the boundary, 0 at all other dofs.

5

I have to check this out. Thank you, sir.

6

It works amazingly nice. Finally I could operate with normal derivative in the same way as any other function.

def normal_for_square_mesh():
    mesh = UnitSquareMesh(100, 100)
    n = FacetNormal(mesh)
    V = VectorFunctionSpace(mesh, "CG", 1)
    u = TrialFunction(V)
    v = TestFunction(V)
    a = inner(u, v) * ds
    l = inner(n, v) * ds
    A = assemble(a, keep_diagonal=True)
    L = assemble(l)

    A.ident_zeros()
    nh = Function(V)

    solve(A, nh.vector(), L)
    File("nh.xml") << nh
    return nh


normal_function = []


def get_normal_derivative(function):
    if not normal_function:
        normal_function.append(normal_for_square_mesh())
    normal = normal_function[0]
    mesh = UnitSquareMesh(100, 100)
    f = FunctionSpace(mesh, "CG", 1)
    function = interpolate(function, f)
    V = VectorFunctionSpace(mesh, "CG", 1)
    return project(inner(normal, project(grad(project(function, f)), V)), f)