Defining a function from data

Question

I am looking to define a Robin boundary condition for the (screened) Poisson equation:

$$
-\Delta u + \kappa^2u = f \text{ in } \Omega
$$$$
\beta u + \frac{\partial u}{\partial n} = 0 \text{ in } \partial \Omega,
$$

where $\beta$ is a function defined on the boundary. For my given mesh, I know the values
of $\beta$ at the nodes. Can I define the Robin BC directly, without interpolating $\beta$? If not, any reference to documentation interpolation on the boundary will be most helpful.
In the code I'm using, I just use $\beta \equiv 1$. This is a (slight) modification of a piece of code written by John Burkardt.

from dolfin import *
import numpy as np

m = n = 500
kappa = 0.3

# I want to use these betas that I retrieve from memory.
# the size of betas is exactly the size of the boundary.
betas = np.load("data/" + filename + ".npy")

# But now I just use this constant beta
beta = 1.0 

mesh = UnitSquareMesh( m, n )

V = FunctionSpace ( mesh, "CG", 1 )
u = TrialFunction ( V )
v = TestFunction ( V )

#  Define the bilinear form, the left hand side of the FEM system.
AUV = inner ( grad ( u ), grad ( v ) ) * dx + kappa*kappa*u*v*dx  + beta*u*v*ds

#  Define the linear form, the right hand side. Later we add point sources.
LV = Constant ( 0.0 ) * v * dx

#  Assemble the system matrix and right hand side vector, 
A, B = assemble_system ( AUV , LV )

#  Modify the right hand side vector to account for point source(s)
delta = PointSource ( V, Point ( 0.05, 0.5 ), 10.0 )
delta.apply ( B )

#  Solve the linear system for u.
u = Function ( V )
solve ( A, u.vector(), B )

Comments

Just a suggestion, but since you know the values at the boundary nodes, can't you define a Function B and then do something like B.vector()[:] = known_node_values. Here known_node_values can be zero wherever a node is not on the boundary, since it doesn't matter anyway.

---

@maartent Thank you for your comment. Following your advice, I started looking in this direction. However, wouldn't it be more natural to define a function that lives just on the boundary? Any further references would be appreciated.

---

I agree it would be more natural, but as far as I know using functions defined on the boundary in your boundary conditions is something yet to be implemented. You can search for 'Submesh' and/or 'Boundarymesh' to find other questions that might be helpful to you (or not).

---

So you suggest doing something like:

known_node_values = get_node_values_somehow()   

mesh = UnitSquareMesh( 50, 50 )

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

B = Function( V )

B.vector()[:] = known_node_values

Any tips on how to align my known values with the nodes?
Shouldn't the last line be like this?

B.vector().array()[:] = known_node_values_as_numpy_array

---

B.vector()[:] = rhs should work, when rhs is a numpy array of the correct size.

Where did you get these node values from? From another problem that you solved in FEniCS (in which case there is probably a way to map them to your new function -- see previous comment) or from a table or so (in which case you can eg. subclass an Expression B instead of a Function, and define it such that it returns the correct values at the boundary)?

---

I got it from some optimization problem I solve for beta in python. My problem is that the correspondence between what I would consider a natural order for the unit square does not agree with FEniCS's ordering.

Also - how would I go about using an Expression? Transform the numpy array to a text file and add some C code to make sense of it?

---

http://fenicsproject.org/documentation/dolfin/dev/python/programmers-reference/functions/expression/Expression.html

should help you on your way (you don't necessarily need C code). You can add an __init__ method where you load your numpy array, and in your eval method you can then use this array, perhaps in combination with an interpolator from scipy.
This has the advantage that is is more general. But since you claim you have the values for the nodes on the boundary already (in combination with their coordinates, I assume), it might be worthwhile to find the corresponding nodenumbers with mesh.coordinates(), and then manually applying the bc's by altering your matrix and rhs. I don't have to much experience with the latter though, so you should check the documentation for this.