I have a diffusion problem which can be broken down to be:

$-\Delta u = f(u) $ on $\Omega ~/~ \Omega_{int}$

$u = 1$ on $\Omega_{int}$

Note that this is an internal Dirichlet constraint to the value 1.0 in the center circle, unlike the standard textbook examples with homogeneous boundary conditions. (see picture) Think of it like a source in the middle, which diffuses into the domain with f(u) = u(1-u) being a logistic growthterm. I want that reaction to stop at a certain maximum total mass.

I want to impose an integral condition which will stabilize the solution:

$P(u) = \int_{\Omega} u ~ dx =V_{tot}$

Now, I tried to read into the problem, and I find that it is possible to add penalty terms to the bilinear form of the finite element discretization, but I am deeply puzzled as to how to properly define the penalty term so that my Newton solver will behave correctly.

Multiplying with a test function v, integrating over the domain we may formulate the problem without the volume constraint like this:

find u, v so that

$ r(u,v) = a(u,v)+ l(v) = 0 $

where $r(u,v)$ is a global residuum (To be minimized by the Newton solver), $a(u,v)$ is the canonical discretisation of the diffusive term and $l(v)$ the disc. of the reaction term.

Now in the language of Lagrange multipliers, we have a functional $r(u,v)$ and a constraint $P(u)$ on the volume so that the new Functional would look like:

$L(u,v,\lambda)= a(u,v)+ l(v) -\lambda(P(u)-V_{tot}) = 0 $

  • How do I apply this to my Finite Element Method in a clean way?
  • How would I determine $\lambda$?
  • Is there any readable literature on how to construct volume-constraints via penalty terms?
  • Is there a more readable/hands-on literature/Paper on Lagrange multipliers in the context of FE-Methods than the work of Ivo Babuska?
  • How do I make sure that the penalty term, which will probably contain a domain-wide integral will not mess up my Newtonian (since the necessary calculations for the Jacobian are no longer local)?

enter image description here

  • 1
    $\begingroup$ If you have $u=1$ on $\Omega_{int}$, then the solution is equivalent to the one of the Laplace equation on $\Omega-\bar\Omega_{int}$ with boundary conditions equal to 1 on $\partial\Omega_{int}$. This might be the easier way to solve the problem -- just cut out $\Omega_{int}$. $\endgroup$ Mar 18 '19 at 17:45
  • $\begingroup$ Hi. I don't think that's it. If I do not constrain the total volume, then the domain would simply 'fill up' diffusively until the solution reaches u=1 everywhere. Note that I keep the outer boundary of the box unconstrained! Thanks for answering! $\endgroup$
    – MPIchael
    Mar 19 '19 at 8:45
  • $\begingroup$ Well, if $u=1$ in $\Omega_{int}$, then $u=1$ on $\partial \Omega_{int}$. And if $-\Delta u=f$ in $\Omega \cup \Omega_{int}$, then also $-\Delta u=f$ in $\Omega - \Omega_{int}$. $\endgroup$ Mar 19 '19 at 13:53
  • $\begingroup$ I'll also note that you can't have $-\Delta u=f$ in $\Omega_{int}$ and $u=1$ in $\Omega_{int}$ at the same time, unless $f=0$ in $\Omega_{int}$. $\endgroup$ Mar 19 '19 at 13:54
  • 1
    $\begingroup$ @MPIchael assuming you intend your solution to be unique, it's likely you do implicitly have a boundary condition on the outer domain, with the most likely candidate being a "no-flux" condition like $\mathbf{n}\cdot\nabla u = 0$. Otherwise there are multiple solutions with a given total volume. $\endgroup$
    – origimbo
    Mar 19 '19 at 19:09

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.