dolfin team mailing list archive

[Anders Logg <logg@xxxxxxxxx>]

On Mon, Feb 11, 2008 at 04:02:29PM -0600, Matthew Knepley wrote:
> On Feb 11, 2008 1:57 PM, Matthew Knepley <knepley@xxxxxxxxx> wrote:
> > On Feb 11, 2008 10:27 AM, Anders Logg <logg@xxxxxxxxx> wrote:
> > >
> > > On Mon, Feb 11, 2008 at 08:05:57AM -0600, Matthew Knepley wrote:
> > > > On Feb 11, 2008 5:05 AM, Anders Logg <logg@xxxxxxxxx> wrote:
> > > > > On Mon, Feb 11, 2008 at 11:54:14AM +0100, Kristen Kaasbjerg wrote:
> > > > > > Hi dolfin users
> > > > > >
> > > > > > I am trying to use dolfin for the following simple
> > > > > > electrostatic problem:
> > > > > > - finding the electrostatic potential in a domain composed
> > > > > >   of two subdomains with different dielectric constants.
> > > > > >   This amounts to solving Laplace equation with appropriate
> > > > > >    BC's.
> > > > > >
> > > > > > One thing that I cannot figure out how to do is to tell dolfin
> > > > > > that I have 2 subdomains and on the internal boundary between
> > > > > > these two domains, the normal derivative of the potential has
> > > > > > a discontinuity given by the difference between the dielectric
> > > > > > constants.
> > > > > > Could anyone give me a hint of how to approach this problem
> > > > > > with dolfin ? I am relatively unexperienced to FEM so a nice
> > > > > > reference would also be welcome.
> > > > > >
> > > > > > Regards
> > > > > > Kristen
> > > > >
> > > > > The easiest thing to do (if there is just a coefficient in your
> > > > > problem that is discontinuous) is to just define a Function for the
> > > > > dielectric constant and make it discontinuous. Then just plug it in to
> > > > > your equation.
> > > > >
> > > > > Look in the demos (under src/demo/pde) for how to define a function.
> > > > >
> > > > > For example, if some parameter is 0 or 1 depending on whether x is
> > > > > below or above 0.5, then just do something like this in eval():
> > > > >
> > > > >   if (x[0] > 0.5)
> > > > >     return 1.0;
> > > > >   else
> > > > >     return 0.0;
> > > >
> > > > This is not going to give the correct answer unless you do the boundary integral
> > > > that comes from the integration by parts. Or am I missing something?
> > > >
> > > >    Matt
> > >
> > > I don't know. Maybe I'm ignorant but if you have -div(a*grad(u)) = f,
> > > then it looks to me that you may just go ahead and integrate by parts
> > > even if the coefficient a is discontinuous. grad(u) is already
> > > discontinuous (for the discretization) so I don't see the difference.
> >
> > I do not agree. There is a jump term along the boundary. You can get
> > some solution without it, which looks alright, it is just wrong.
> 
> How about this example. You can tell me what is wrong:
> 
>   Take a 1-D domain from -1 to 1 which is divided at 0 into
> \epsilon_0 and \epsilon_1. Then we have something like
> 
>   \int^1_{-1} v div (\epsilon grad u)
> 
>  = \int^0_{-1} v div (\epsilon_0 grad u) + \int^1_0 v div (\epsilon_1 grad u)
> 
>  = -\int^0_{-1} grad v \epsilon_0 grad u + \epsilon_0 v grad u |_0
>      -\int^1_0 grad v \epsilon_1 grad u - \epsilon_1 v grad u |_0
> 
>  = -\int^1_{-1} \epsilon grad v grad u - (\epsilon_1 - \epsilon_0) v grad u |_0
> 
> Dmitry is now of the opinion that this is wrong based on distribution theory
> arguments. However, I was sure that physically we get a layer of induced
> charge at a dielectric jump proportional to the jump times the field. I will try
> and look it up in Jackson.
> 
>    Matt

I would say that what you want to solve is

    min 1/2 a(v, v) - (v, f)

where a(v, v) = \int epsilon grad v . grad v. This is well-defined for
discontinuous epsilon. The solution of the minimization problem is
just the function u that solves

   a(v, u) = (v, f) for all v

which is also well-defined for discontinuous epsilon.

If you go ahead and take div of epsilon grad(u) then you need to
interpret div of something discontinuous appropriately, and I suspect
the appropriate way to interpret it is to add a jump term that will
exactly cancel out those terms that you get when integrating by parts.

-- 
Anders