Subject: [Fenics] 3D Diffusion Problem with Fenics (Python Interface)

From: Lorenzo Isella <[email-masked]>
Date: Thu, 03 Mar 2011 11:55:59 +0100

Dear All,
I actually posted something along this lines quite some time ago.
I apologize in advance for the long, detailed message I am posting, but 
any help about any step involved, would help me a lot.
What I am looking for should be, I am sure, something rather simple for 
somebody knowledgeable about Fenics. It is a fairly general problem, but 
I will try to make it more concrete with an example (alas, no code, 
otherwise I would have already solved this myself).
Consider a cube with size L=50 and centered about the origin of your 
coordinate system (i.e. every side stretches along the interval [-25,25] 
along the x,y, and z axis).
Then consider a simple straight linear chain consisting of three spheres 
having radius R=1. The y and z coordinates of the (centers of) the 
spheres are precisely y=z=0) whereas their x coordinates are x_1=-2, 
x_2=0 and x_3=2.
To make a long story short, I have a straight linear chain made up of 3 
touching impenetrable spheres oriented along the x axis and with 
center-of-mass in the origin of my coordinate system (0,0,0).
Now, the equation I would like to solve is nothing terrible, but simply

D\nabla^2\rho=0 (where \rho=\rho(x,y,z) is scalar you can think of a a 
density of some kind), with boundary conditions

rho=1 along the facets of the cube and rho=0 along the surface of the 
spheres (I can set D=1 if I want to). Nothing needs to be solved within 
the 3 spheres. I can see the problem in at least two equivalent ways
(1) the computational domain is given by the cube "minus" the three 
spheres, i.e. it is as if I was carving the spheres inside the cube or
(2) I actually have D=0 inside the 3 spheres (hence nothing diffuses 
inside them) with the boundary condition \rho=0 on their surface.

In any case, there is nothing I need to compute within the spheres.

The final goal is to be able to integrate the diffusional flux (or its 
absolute value) along the surfaces of all the spheres, so I just need a 
reliable number out of all of this. If possible, I would like to do 
everything with a single Python script, without resorting to an external 
tool to generate the mesh.

If I manage to get this example up and running in Fenics, then I am 
almost done (i.e. it is conceptually nothing difficult to enlarge the 
cube, separate the spheres or add more spheres), but I also need some 
way to tune the mesh explicitly to ensure that my results are mesh 
independent.
Any help from the list is really appreciated.
Many thanks

Lorenzo