Question on Dolfin poisson equation for source measurement on borders

Legacy DOLFIN
,
1

Hello,

I want to know if it is possible, using the solution (u) with a source term (f) of the poisson equation, to measure that solution on the boundaries instead of an “inner” measure.

It would be similar to the Electrical Impedance Tomography measure, where we have eletroces on the boundary measuring the solution (u) :

image

Here is my code inspired mainly from the dolfin example Pure neumann poisson :

from dolfin import *
from mshr import *

resolution = 15
domain_geometry = Circle(dolfin.Point(0.5, 0.5), 0.5+0.01)
mesh = generate_mesh(domain_geometry, resolution)

# Build function space with Lagrange multiplier
P1 = FiniteElement("Lagrange", mesh.ufl_cell(), 1)
R = FiniteElement("Real", mesh.ufl_cell(), 0)
W = FunctionSpace(mesh, P1 * R)

# Define variational problem
(u, c) = TrialFunction(W)
(v, d) = TestFunctions(W)
f = Expression("10*exp(-(pow(x[0] - 0.5, 2) + pow(x[1] - 0.5, 2)) / 0.02)", degree=2)
g = Constant("0.0")
a = (inner(grad(u), grad(v)) + c*v + u*d)*dx
L = f*v*dx + g*v*ds

# Compute solution
w = Function(W)
solve(a == L, w)
(u, c) = w.split()

# Save solution in VTK format
file = File("neumann_poisson.pvd")
file << u

Another code I did based on this poisson example Poisson equation :

from dolfin import *
from mshr import *

# Create mesh and define function space
resolution = 15
domain_geometry = Circle(dolfin.Point(0.5, 0.5), 0.5+0.01)
mesh = generate_mesh(domain_geometry, resolution)
V = FunctionSpace(mesh, "Lagrange", 1)

def boundary(x, on_boundary):
    return on_boundary

# Define boundary condition
u0 = Constant(0.0)
bc = DirichletBC(V, u0, boundary)

# Define variational problem
u = TrialFunction(V)
v = TestFunction(V)
f = Expression("10*exp(-(pow(x[0] - 0.5, 2) + pow(x[1] - 0.5, 2)) / 0.02)", degree=2)
g = Constant(0.0)
a = inner(grad(u), grad(v))*dx
L = f*v*dx + g*v*ds

# Compute solution
u = Function(V)
solve(a == L, u, bc)

# Save solution in VTK format
file = File("poisson.pvd")
file << u

Thank you for your help.

Great day to you :slight_smile:

4

Unfortunately I don’t quite understand your question. Could you please try to be more concise, and maybe add the mathematical formulation of what you want to achieve.

To me it seems like you want to compute (u, f)_{l^2(\mathbb{R}^N)} for only nodes on the boundary?
or do you want to compute (u, f)_{L^2(\partial\Omega)} ?

5

Hello Dokken,

Thank you for your reply.

A first approach would be to compute the 2nd formulation you gave, which means that my measure will be continuous on the entire boundary.

After that, I would like to achieve your 1st formulation which is (for example) specifying number of nodes to be equivalent to one electrode, and do a step jump for the other electrode, in order to compute that 1st formulation.

Thank you and have a great day.

6

I dont see the problem of computing these.

Am I missing something?

For 2 you could just use dot(u,f)*ds assuming u and f are dolfin/ufl objects.

For 1 you would extract the coordinates you would like to use (if they align with dof coordinates things become easier), and evaluate u(x_i)f(x_i).

7

I can try to figure this out on my own, but I am not sure if it will be simple.

Do you have an example on how to achieve this ? Or maybe a general approach for both formulations ?

Thank you for your time.

8

Hello again Dokken,

I am not sure that I can adapt your two formulations in a code.

Can you please give me a hint on where to look in order to modify one of the previous codes ?

Thank you.

9

For the L2 formulation on the boundary, just use assemble(inner(u,f)*ds(marked_boundary)).

Please note that this is an open-source, volunteer forum, and there are many posts here every day. I simply do not have the bandwidth to reply to all of them.

For the l2 formulation, it requires some more work, as you would first have to identify what degrees of freedom that are located on your boundary of interest, then extract those from u and f and multiply the arrays.

11

You should send in a tuple of integers that is used in marked boundary to mark the relevant portions, ie if marked boundary contains the markers (1,2,3), and you want to integrate over 1 and 3 you send in (1, 3)

12

Sorry but I can’t figure out what you are saying. I have only one domain wich is my circle mesh, and I want to to solve “u” only on the boundary. How can multiple subdomains be present ?

13

A MeshFunction object can contain more than a single value splitting the facets into sets

As you want to use the whole external boundary, you can evaluate two known function u, F on that boundary, you can simply use assemble(inner(u,f)*ds)

You are using the Word “solve” when I think you mean evaluate?

14

Okey I see, so there is no need for these lines :

########## Added portion ###########################
marked_boundary = MeshFunction("size_t", mesh, mesh.topology().dim()-1)
marked_boundary.set_all(0)

ds = Measure('ds', domain=mesh, subdomain_data=marked_boundary)
####################################################

To reply to your question : For me the word “solve” is related to the fact that I am computing a function u, with consideration to a boundary condition, using a weak formulation (i.e : a == L). Isn’t it the same if we say that we are evaluating that u function over a specific domain (which is our mesh) ?

And I tried your solution, keeping the same code and removing those lines I mentionned above, with :

L = assemble(inner(u, f)*ds)

I get the following message followed by error :

No Jacobian form specified for nonlinear variational problem.
Differentiating residual form F to obtain Jacobian J = F.

And the error message specifies :

Error : Unable to define nonlinear variational problem F(u; v) = 0 for all v.
Reason : Expecting the residual F to be a linear form (not rank 2).

Thank you.

15

I have no clue what you are trying to do any longer.

What i thought you wanted to:

  1. solve a linear problem a(u,v)=L(v) for all v in some function space
  2. Evaluate the solution u (after solving) at the boundary with the inner product defined on the boundary
    Could you please write out the mathematical formulation of what you are trying to do if the above steps are not what you are trying to do.
16

My bad, I mixed the 2 terms. That is exactly what I want to achieve.

17

If we impose g = Constant(0.0) in the variational problem, then we will always measure assemble(inner(u,f)*ds) equal to 0, no ?

If not, then I don’t know how I can evaluate the solution u (after solving) at the boundary.

Thank you

18

Why?
You are emposing du/dn=g on the boundary, which does not make that integral (ufds) zero.

I don understand your reasoning here.

A finite element function can be evaluated anywhere in the domain, as it is defined as u(x)=\sum_{i=0}^N c_i \phi_i(x), where c_i are the coefficients/dofs obtained by calling solve.

1 Like
20

You set a Dirichlet BC on the outer boundary, hence the solution u is zero and its boundary integral (multiplication by f doesn’t really matter) will be zero.

1 Like
22

You are asking to evaluate the boundary integral of f times u, where both f and u are scalar fields, so the result is a (single) scalar.

1 Like
23

I understand, so how can I evaluate u on boundary ? (I should have exact amount of values than u.vector().get_local(), with respect to the source term f)

Edit : I don’t think that the amount of values will be the same, since we only evaluate on the boundary, and u.vector().get_local() gives u values for both boundary and inner. So, the amout of u values on boundary should be less.

24

A simple way could be to use V.tabulate_dof_coordinates() to get the index of a DOF and its associated coordinate. Once you have all coordinates, you can filter out those which are on the boundary and those which are in the interior.

There may be more elaborate ways of doing that, but I won’t have the time to prepare an example. Those ways may have already been discussed in this forum, search e.g. using MeshView (or even multiphenics, which I used to develop).

25

I tried a method with a dummy example. i.e :

from dolfin import *
mesh = UnitSquareMesh(4, 4)

V = FunctionSpace(mesh, 'CG', 1)
bc = DirichletBC(V, 1, "on_boundary")
u = Function(V)
bc.apply(u.vector())

plot(mesh)

In order to obtain the boundary coordinates, I just do :

bmesh = BoundaryMesh(mesh, "exterior", True)
print(bmesh.coordinates())

Now the goal is to evaluate the u solution on those coordinates.

I wanted to make sure that if I print the value of the u solution on specific coordinates, then i’ll get the value I have in print(u.vector().get_local()) :

image

As you can see, starting from 0 index, the 4th u value is equal to 0.

Now if we get the 4th coordinates (starting from 0 index) through print(mesh.coordinates()) :

image

We can see that print(u(1, 0)) is equal to 1 and not 0 !

image

Is this related to this post Dof and coordinates, where MiroK specified that : there is a difference between how coordinates of vertices in mesh are ordered and the order of degrees of freedom of functions from functions spaces over that mesh ?

Of course, my question is still related to the fact that I want to evaluate u on boundary coordinates that I generated (so I would get like an array of u values for those coordinates).

Thank you for your time.