How to apply a robin boundary condition using a solved variational problem

Question

Hi all,

In my current code I first solve the poisson equation in a large mesh with different conducting properties.
Next I want to solve in a submesh which is much smaller and will be changed to contain microscopic properties.

I would like to be able to find robin boundary conditions that I can apply on boundaries of the submesh based on the original calculation.

β u + ∂u / ∂n = 0 in ∂Ω

So I have to find a β that can "simulate" the rest of the domain. If I keep the submesh unchanged it should result in the same solution by applying the same flux boundary condition as we find in the "total" solution

How could I compute these β values?

Answer 1

Yes this is possible. You can interpolate the results of your FE solution on the initial mesh to be used in the boundary computation of the mesh of the subdomain at a geometric location.

For example, you can interpolate the value of a function at the point (0.5, 0.5) as follows:

expr = Expression('sin(pi*x[0])*sin(pi*x[1])', element=V.ufl_element())
u = Function(V)
u.interpolate(expr)
print u(0.5, 0.5)

Comments

Thank you for answering.
But my problem description might not be clear enough, or I don't understand how to apply your solution to my problem. Here is some code where I describe in more detail at what point I am having difficulty.

from dolfin import *

''' Create mesh and define function space'''

mesh = UnitSquareMesh(5, 5)

V = FunctionSpace(mesh, 'CG', 1)

''' Define Dirichlet boundary conditions'''

u0 = Expression('1 + x[0]x[0] + 2x[1]*x[1]', degree=2)

class DirichletBoundary(SubDomain):

    def inside(self, x, on_boundary):

        tol = 1E-14   # tolerance for coordinate comparisons

        return on_boundary and \

               (abs(x[0]) < tol or abs(x[0] - 1) < tol)

u0_boundary = DirichletBoundary()

bc = DirichletBC(V, u0, u0_boundary)

''' Define variational problem'''

u = TrialFunction(V)

v = TestFunction(V)

f = Constant(-6.0)

g = Expression('-4x[1]', degree=1)

a = inner(grad(u), grad(v))dx

L = fvdx - gvds

''' Compute solution'''

sol = Function(V)

solve(a == L, sol, bc)

This first part is just taken from the tutorial. Next I create a submesh and want to find robin boundaries that simulate the rest of the mesh:

''' create submesh and functionspace'''

class MidSquare(SubDomain):

    def inside(self, x, on_boundary):

        return (between(x[0], (0.4, 0.6)) and between(x[1], (0.4, 0.6)))

subdomain = CellFunction('size_t', mesh)

subdomain.set_all(0)

midsquare = MidSquare()

midsquare.mark(subdomain,1)

submesh = SubMesh(mesh, subdomain, 1)

subV = FunctionSpace(submesh, 'CG', 1)

dx_mydomains = CellFunction('size_t', submesh)

dx_mydomains.set_all(0)

ds_myfacets = FacetFunction('size_t', submesh)

ds_myfacets.set_all(0)

'''Define variational problem'''

u = TrialFunction(subV)

v = TestFunction(subV)

f = Constant(-6.0)

'''

Do something to figure out the correct betas here based on the first solution

It should be something like n.grad(sol) = beta * u along the boundary between mesh and submesh

The beta values should make sure we find the same solution in the submesh,

without having to solve the whole mesh

I tried something like below, but how should i apply that as a boundary condition?

n=FacetNormal(mesh)

dvdn=dot(nabla_grad(sol),n)

beta=dvdn/sol

'''

beta = Expression('x[0]', degree=1) #example beta

a = inner(grad(u), grad(v))dx(subdomain_data=dx_mydomains) + betauvds(subdomain_data=ds_myfacets)

L = fvdx(subdomain_data=dx_mydomains)

'''Compute solution'''

subsol = Function(subV)

A, b = assemble_system(a, L )

solve(A, subsol.vector(), b)