# Poisson equation the stubborn way

**URL:** <https://fenicsproject.discourse.group/t/poisson-equation-the-stubborn-way/1645>\
**Category:** variational formulation\
**Created:** [September 30, 2019, 7:17pm UTC](https://fenicsproject.discourse.group/t/poisson-equation-the-stubborn-way/1645 "2019-09-30T19:17:49Z")\
**Posts on this page:** 5\
**Page:** 1

<div class="post-metadata">

**Author:** ![jw3126](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/jw3126/32/508_2.png) [@jw3126](https://fenicsproject.discourse.group/u/jw3126)\
**Post date:** [September 30, 2019, 7:17pm UTC](https://fenicsproject.discourse.group/t/poisson-equation-the-stubborn-way/1645/1 "2019-09-30T19:17:49Z")

</div>

I am trying to understand what is a good variational formulation and what is a bad variational formulation and why. It seems common knowledge, that one should minimize the degree of the highest differential operator. Here is an example of a naive Poisson equation where I did not do that:

```python
from fenics import *

mesh = IntervalMesh(10, -1.0, 1.0)
V = FunctionSpace(mesh, "P", 1)
phi = Function(V)
rho = Constant(1.0)
v = TestFunction(V)
F = inner(-div(grad(phi)), v) * dx - inner(rho, v) * dx # does not work
# F = inner(-grad(phi), -grad(v)) * dx - inner(rho, v) * dx # works
bc = DirichletBC(V, Constant(0.0), "on_boundary")
solve(F == 0, phi, bc)

```

It throws an error:

```auto
*** -------------------------------------------------------------------------
*** Error: Unable to solve nonlinear system with NewtonSolver.
*** Reason: Newton solver did not converge because maximum number of iterations reached.
*** Where: This error was encountered inside NewtonSolver.cpp.
*** Process: 0
*** 
*** DOLFIN version: 2019.1.0
*** Git changeset: unknown
*** -------------------------------------------------------------------------

```

It even throws an error, if I give it an analytic solution as initial guess:

```python
phi = interpolate(Expression("-0.5*(x[0]-1)*(x[0]+1)", degree=2), V)

```

- Why this error?
- Why is `inner(grad(phi), grad(v))` better then `inner(-div(grad(phi)), v)`?
- Can I change the definition of the vector space `V` to make `inner(-div(grad(phi)), v)` work?

---

<div class="post-metadata">

**Author:** ![nate](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/nate/32/17_2.png) [@nate](https://fenicsproject.discourse.group/u/nate)\
**Post date:** [September 30, 2019, 8:05pm UTC](https://fenicsproject.discourse.group/t/poisson-equation-the-stubborn-way/1645/2 "2019-09-30T20:05:31Z")

</div>

Someone else may answer your question specifically, but for a grasp of the fundamentals I recommend reading the introductory literature.

One of many which I like:  
[The Finite Element Method: Linear Static and Dynamic Finite Element Analysis By Thomas J. R. Hughes](https://books.google.com/books?id=cHH2n_qBK0IC&printsec=frontcover&source=gbs_ge_summary_r&cad=0#v=onepage&q&f=false)

Another I like:  
Numerical Solution of Partial Differential Equations by the Finite Element Method by Claes Johnson

---

<div class="post-metadata">

**Author:** ![kamensky](https://avatars.discourse-cdn.com/v4/letter/k/e95f7d/32.png) [@kamensky](https://fenicsproject.discourse.group/u/kamensky)\
**Post date:** [October 1, 2019, 2:19pm UTC](https://fenicsproject.discourse.group/t/poisson-equation-the-stubborn-way/1645/3 "2019-10-01T14:19:20Z")

</div>

I second Nate’s suggestion to take some time with the fundamentals. To answer the questions specifically,

- _Why this error?_ Because `dx` is really a sum of integrals over interiors of elements. You can see that, on the interior of each element, -\nabla\cdot(\nabla\phi) is zero when \phi is piecewise linear. Thus, the left-hand side matrix is singular.

- _Why is `inner(grad(phi), grad(v))` better?_ Because the other definition is not correct! Glossing over some technicalities about the theory of distributions, recall that the derivative of a Heaviside function is a Dirac delta. Considering `-div(grad(phi))` on the element interiors only (as per the `dx` measure) is missing the Dirac deltas that should be part of the derivative of the discontinuous function `grad(phi)`. See the following example, where `dS` is a sum of integrals over interior mesh facets:

```python
from fenics import *

mesh = IntervalMesh(10, -1.0, 1.0)
V = FunctionSpace(mesh, "P", 1)
phi = Function(V)
rho = Constant(1.0)
v = TestFunction(V)
F = inner(-div(grad(phi)), v) * dx - inner(rho, v) * dx # does not work

# Because grad(phi) is discontinuous, its derivative includes Dirac deltas
# which must be explicitly added, because *dx is a sum of integrals
# over element interiors.
n = FacetNormal(mesh)
F += inner(jump(grad(phi),n),avg(v))*dS

bc = DirichletBC(V, Constant(0.0), "on_boundary")
solve(F == 0, phi, bc)

# Compare with standard formulation:
phi2 = Function(V)
F = inner(-grad(phi2), -grad(v)) * dx - inner(rho, v) * dx # works
bc = DirichletBC(V, Constant(0.0), "on_boundary")
solve(F == 0, phi2, bc)

# Solutions are identical:
phiDiff = phi - phi2
print(assemble((phiDiff**2)*dx))

```

- _Can I change the definition of the vector space `V` to make `inner(-div(grad(phi)), v)` work?_ Mathematically, yes, you could use a C^1 space of functions, in which case the jump terms in the above example vanish, although this is not supported directly in FEniCS. (You can do it using a [library](https://github.com/david-kamensky/tIGAr) I worked on, though.) However, it would then be exactly equivalent to `inner(grad(phi),grad(v))` (with Dirichlet BCs).

---

<div class="post-metadata">

**Author:** ![jw3126](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/jw3126/32/508_2.png) [@jw3126](https://fenicsproject.discourse.group/u/jw3126)\
**Post date:** [October 2, 2019, 4:37pm UTC](https://fenicsproject.discourse.group/t/poisson-equation-the-stubborn-way/1645/4 "2019-10-02T16:37:37Z")

</div>

Thanks this answer was quite enlightening. Is the fact that `dx` only integrates of the interiors of elements documented somewhere?

---

<div class="post-metadata">

**Author:** ![kamensky](https://avatars.discourse-cdn.com/v4/letter/k/e95f7d/32.png) [@kamensky](https://fenicsproject.discourse.group/u/kamensky)\
**Post date:** [October 2, 2019, 5:23pm UTC](https://fenicsproject.discourse.group/t/poisson-equation-the-stubborn-way/1645/5 "2019-10-02T17:23:22Z")

</div>

It comes down to terminology, since “integrates” is typically understood in the sense of Lebesgue when formulating finite element methods, so integrating over element interiors is equivalent to integrating over the full domain, and it’s simply invalid to talk about integrating an object like -\nabla\cdot(\nabla\phi^h), since it’s not a function. (It’s arguably an abuse of notation to write things like \int \delta\_0(x)v(x)\,dx = v(0); really it’s the application of a linear functional \delta\_0 to the function v, i.e., \langle \delta\_0, v\rangle.) From that perspective, it’s a user error to ask for an ill-defined integral. However, that’s not the full story, since many residual-based stabilized methods actually do include terms that are only integrable on element interiors and should technically be written using a sum over element integrals, which is one reason I can think of why things like `inner(-div(grad(phi)),v)*dx` are allowed to compile without errors.

EDIT: To answer the question more directly, yes, the [UFL manual](https://buildmedia.readthedocs.org/media/pdf/fenics-ufl/latest/fenics-ufl.pdf) explicitly defines the `dx` measure as a sum over element integrals. See the formula at the bottom of page 4.