Solve PDEs without transforming to variational formulation

variational formulation
1

Dear Community members!

I started to learn FEniCS using the book by Langtangen and Logg.
When solving the Poisson-equation, it is suggested, to reformulate the equation

\Delta u = -f

to the variational formulation, i.e. multiply by v \in H^1_0 and integrate by parts to get

\int \nabla u \cdot \nabla v \text{dx} = \int f v \text{dx}.

Just to play around with the code, I wanted to implement this equation instead, that one obtains by not doing the integration by parts:

\int \Delta u \cdot v \text{dx} = -\int f v \text{dx}.

In my understanding u \in H^1 satisfying this expression would still be a weak solution.

The code looks like this:

# spaces and mesh
mesh = UnitSquareMesh(8, 8)
V = FunctionSpace(mesh, 'Lagrange', 1)
u = TrialFunction(V)
v = TestFunction(V)

# BC
u_D = Expression('1', degree=0)
def boundary(x, on_boundary):
    return on_boundary

bc = DirichletBC(V, u_D, boundary)

# equations
f = Constant('-10')
a = dot(div(grad(u)), v)*dx
L = f*v*dx
u = Function(V)
solve (a==L, u, bc)

# plot
plot(mesh)
plot(u)

And produces the following error message:

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-67-38954b2d310f> in <module>
     21 # plot
     22 plot(mesh)
---> 23 plot(u)

/usr/local/lib/python3.6/dist-packages/dolfin/common/plotting.py in plot(object, *args, **kwargs)
    436     # Plot
    437     if backend == "matplotlib":
--> 438         return _plot_matplotlib(object, mesh, kwargs)
    439     elif backend == "x3dom":
    440         return _plot_x3dom(object, kwargs)

/usr/local/lib/python3.6/dist-packages/dolfin/common/plotting.py in _plot_matplotlib(obj, mesh, kwargs)
    302 
    303     if isinstance(obj, cpp.function.Function):
--> 304         return mplot_function(ax, obj, **kwargs)
    305     elif isinstance(obj, cpp.function.Expression):
    306         return mplot_expression(ax, obj, mesh, **kwargs)

/usr/local/lib/python3.6/dist-packages/dolfin/common/plotting.py in mplot_function(ax, f, **kwargs)
    150             if mode == "contourf":
    151                 levels = kwargs.pop("levels", 40)
--> 152                 return ax.tricontourf(mesh2triang(mesh), C, levels, **kwargs)
    153             elif mode == "color":
    154                 shading = kwargs.pop("shading", "gouraud")

/usr/local/lib/python3.6/dist-packages/matplotlib/tri/tricontour.py in tricontourf(ax, *args, **kwargs)
    273 def tricontourf(ax, *args, **kwargs):
    274     kwargs['filled'] = True
--> 275     return TriContourSet(ax, *args, **kwargs)
    276 
    277 

/usr/local/lib/python3.6/dist-packages/matplotlib/tri/tricontour.py in __init__(self, ax, *args, **kwargs)
     33         are described in the docstring of `tricontour`.
     34         """
---> 35         ContourSet.__init__(self, ax, *args, **kwargs)
     36 
     37     def _process_args(self, *args, **kwargs):

/usr/local/lib/python3.6/dist-packages/matplotlib/contour.py in __init__(self, ax, levels, filled, linewidths, linestyles, alpha, origin, extent, cmap, colors, norm, vmin, vmax, extend, antialiased, *args, **kwargs)
    853         self._transform = kwargs.pop('transform', None)
    854 
--> 855         kwargs = self._process_args(*args, **kwargs)
    856         self._process_levels()
    857 

/usr/local/lib/python3.6/dist-packages/matplotlib/tri/tricontour.py in _process_args(self, *args, **kwargs)
     45         else:
     46             from matplotlib import _tri
---> 47             tri, z = self._contour_args(args, kwargs)
     48             C = _tri.TriContourGenerator(tri.get_cpp_triangulation(), z)
     49             self._mins = [tri.x.min(), tri.y.min()]

/usr/local/lib/python3.6/dist-packages/matplotlib/tri/tricontour.py in _contour_args(self, args, kwargs)
     92                              'triangulation')
     93         if not np.isfinite(z_check).all():
---> 94             raise ValueError('z array must not contain non-finite values '
     95                              'within the triangulation')
     96 

ValueError: z array must not contain non-finite values within the triangulation

So apparently, the solution has some non-finite values. But, when implementing the left bilinear form as \nabla u \cdot \nabla v the code does produce a solution. I dont understand that behaviour, as the two formulations should be equivalent. Is the problem here, that u \notin \mathbf{C}^2 is somehow assumed in the implementation?
I am more than thankful for any kind of help.

2

I’m guessing some (or maybe all) the coefficients of your solution are nan, which is why you get the error you’re seeing.

Degree 1 Lagrange basis functions are continuous between cells, but do not have a continuous derivative. (ie u\in C^1, but u\not\in C^2 as you said). This means that \Delta u is not well defined on the edges of the cells, so to term

\displaystyle\int \Delta u\cdot v\,\mathrm{d}x

cannot be computed correctly.

To see why this is so, you think think about the 1D case. In 1D, the basis functions are hat functions. The derivatives (\nabla u) of these are step functions. The second derivatives (\Delta u) are therefore delta functions. These delta functions should affect the integral (but will not in the FEniCS implementation).

In order to assemble

\displaystyle\int \Delta u\cdot v\,\mathrm{d}x

you will need to use C^2 continuous basis functions. I think these are not supported in FEniCS; and they are not (yet) supported in FEniCSx.

2 Likes
3

Hi Mathew, thank you for your quick answer!

So FEniCS does not provide C^2 basis functions - that implies, any equation has to be written in a formulation, where strictly u\in C^1. Does that mean for an equation like this

\Delta(\Delta u) = -f

one needs to apply integration by parts several times?

By using a second order function space

V = FunctionSpace(mesh, 'Lagrange', 2)

one would get a basis of second order polynomials, like it is described on your website DefElement: Lagrange. If the basis functions are second order, they are \in C^2 locally on each element.
With these second order elements, I still get the same error from the code, though.
Is the problem then caused by lacking regularity of the basis functions on the border/edge of each element?

4

Hello,

Higher order Lagrange elements are still not C^2 continuous, as it is only enforced that their value is the same on mesh entities shared between two cells.

In order to have C^2 continuity, you need to have DOFs that endure that the derivative is continuous between cells. For example, Hermite elements (DefElement: Hermite) are defined using point evaluations and point evaluations of the derivative.

If you’re interested, the reason we can’t yet support Hermite elements in FEniCSx is due to needed a more general pull back/push forward map, as described in this R. Kirby paper.

1 Like
5

In addition to @mscroggs’s points, you could consider weakly enforcing continuity of the first derivative of your FE solution as part of the C^2 requirement. There’s a very old demo which could easily be updated for dolfinx.

There’s also @kamensky’s tIGAr which offers C^k continuity with a spline basis. E.g. the bihamonic demo.

2 Likes
6

I’ll also mention an old post here, where I answered a similar question in some more detail.

3 Likes
7

Thank you very much everybody!