normal flux across the internal interface

[Erik Svensson]

Hello,

I like to find the electrostatic potential, u, in a domain with two materials (solving the Laplace equation in 2D and 3D). On parts of the boundary I have Dirichlet boundary conditions and on other parts I have homogeneous Neumann boundary conditions. I use different dielectric constants, e1 and e2, in part 1 and 2 of the domain, see the fig below. This is all fine up to this point.

Now I would like to have the normal flux across the internal interface, between the materials, to fulfill the jump condition:

n * (e1*grad (u) – e2*grad(u)) = 0,

where n is the unit normal vector pointing out of one of the domains.
 
This is not already imposed by the method, right?

If not, what is a good way of implementing it in deal.ii?

So far I came across a few hurdles, trying it myself.

How do I label internal faces? Boundary ids can only take a value for external faces, on the external boundary (boundary ids on internal faces are reserved to -1, right?).

I am thinking about adding constraints, calculating the constraints in the assembly loop, But then I can not figure out how to add elements to lines in the constraint matrix, using ‘add_entry’ to the ConstraintMatrix will overwrite existing elements at row i and col j, right? I am missing a += operator, adding a number to row i and col j that already exists. Or perhaps I am on a complete wrong track, misunderstanding etc.

Best regards, Erik

[Wolfgang Bangerth]

On 03/24/2016 01:21 AM, Erik Svensson wrote:
>
> Now I would like to have the normal flux across the internal interface,
> between the materials, to fulfill the jump condition:
>
> n * (e1*grad (u) – e2*grad(u)) = 0,
>
> where n is the unit normal vector pointing out of one of the domains.
>
> This is not already imposed by the method, right?

It's not satisfied exactly by the finite element method, but in the weak sense
it is. There is nothing you need to do explicitly to make it work.

Best
W.

--
------------------------------------------------------------------------
Wolfgang Bangerth email: bang...@math.tamu.edu
www: http://www.math.tamu.edu/~bangerth/

[Weixiong Zheng]

Erik,

Constraints would be possibly one way to go. But I will talk about what I did in weak sense.

If you are using continuous finite element in each of the subdomains, it's not clear how you can impose such a condition. However, what I implemented recently in a similar problem with multiple subdomains that required me imposing such a condition, I go to step-46.

Basically, you can choose to use hp classes, which allows you impose different fe for different subdomain and collect all the different fe's in a container pretending to be a ordinary large fe. At the end of the day, the formulation will be like whatever finite element you'd like to use in each subdomain but have a interface condition performing like numerical flux in DG on the interface to mimic your interface condition.

If you can read through the whole step-46, it shall help. And it will be fun to code such a thing out.

Best,

Weixiong

在 2016年3月24日星期四 UTC-4上午2:21:20，Erik Svensson写道：

[Erik Svensson]

Thank you for your answers!

I was actually thinking in terms of continues finite elements. Thought it would be nice to implement the jump in a strong sense. Perhaps this is not possible. Can't see why, if I involve nodes not only on the interface but also in the interior of the domain (I mean the constraints are not spreading all over the domain -nodes in cells with a face on the interface would perhaps be sufficient to use).

In any case case. For the 'regular' implementation Deal.ii nicely deals with the material discontinuity, see the image above.

But the idea was to impose the jump, some how, such that the method could stick to larger elements at the interface without triggering refinements there.

[Wolfgang Bangerth]

> I was actually thinking in terms of continues finite elements. Thought it
> would be nice to implement the jump in a strong sense.

You can't. Imagine a line where the coefficient on one side is e1 and on the
other side e2. Then you know that

n . (e1 grad(u1)) = n . (e2 grad(u2))

where u1, u2 are the solution function evaluated at the two sides of the
interface.

Now imagine that this line cuts through part of the domain where e1=e2, i.e.,
we're considering a line where the coefficient is not actually discontinuous.
The condition then in essence means that the solution's gradient is the same
from the two sides, i.e., the solution is continuous differentiable (in the
space C^1).

Now think what would happen if you tried to enforce something like this
strongly on a discrete solution. For simplicity, take linear triangles. If you
want the gradient on both sides to be the same, that means that the function
on both sides of the interface is simply linear across the interface. Since
this has to be the case at *all* internal interfaces, you'll end up with a
discrete solution that is *globally* just a linear function -- in general a
rather poor "approximation" of the exact solution, and one that doesn't get
better by mesh refinement.

In other words, you really don't want to enforce the flux condition strongly :-)

> In any case case. For the 'regular' implementation Deal.ii nicely deals with
> the material discontinuity, see the image above.
>
> But the idea was to impose the jump, some how, such that the method could
> stick to larger elements at the interface without triggering refinements there.

That's a separate issue. You are using the Kelly error estimator (I suppose)
which looks at the difference of the gradients on both sides of a cell
interface. What you want to do is use a refinement criterion that looks at the
*fluxes* at both sides of a cell interface, i.e., the gradient times the
coefficient. You can achieve this by passing a coefficient to the Kelly error
estimator class.

[Erik Svensson]

Thanks a lot for detailed answer!

Indeed, I'm using the Kelly refinement.

I would need to think more about what I would like to achieve.

/Erik

[Michael Harmon]

I had to do this same problem.

The way to do this is to use a mixed finite element method with Raviart-Thomas elements as done in step 20.  For electrostatics you should have two conditions across the interface:

1.) the jump in potential (phi) is continuous
2.) the jump in the displacement electric field (e_{i} grad(phi)) is continuous (which is what you wrote)

Both interface conditions will automatically be satisfied weakly using the mixed method.  The second constraint is satisfied by the Raviart-Thomas elements itself. 

Is your dielectric constant just a discontinuous step function or spatially variable within each subdomain?

[Erik Svensson]

I apologize (to Michael Harmon and others who follow this thread) for not writing back, promptly! This forum is great and must be used respectfully -I apologize again! I dropped this thread for some time.

Thanks Michael for the suggestion to use Raviat-Thomas elements!

My dielectric constant i just a discontinuous step function across the material interfaces.

/Erik