# Is it possible to solve original or weak formulations of the boundary value problems in DOLFINx?
**URL:** https://fenicsproject.discourse.group/t/is-it-possible-to-solve-original-or-weak-formulations-of-the-boundary-value-problems-in-dolfinx/16282
**Category:** dolfinx
**Created:** [November 12, 2024, 2:18pm UTC](https://fenicsproject.discourse.group/t/is-it-possible-to-solve-original-or-weak-formulations-of-the-boundary-value-problems-in-dolfinx/16282 "2024-11-12T14:18:11Z")
**Posts on this page:** 6
**Page:** 1
### Author:  [@marbor](https://fenicsproject.discourse.group/u/marbor)
#### Post date: [November 12, 2024, 2:18pm UTC](https://fenicsproject.discourse.group/t/is-it-possible-to-solve-original-or-weak-formulations-of-the-boundary-value-problems-in-dolfinx/16282/1 "2024-11-12T14:18:11Z")
For purely Dirichlet boundary value problem:
\begin{eqnarray} -\nabla^2 u = f \ in \ \Omega \nonumber \\ u = \bar{u} \ on \ \Gamma\_D \nonumber \\ \end{eqnarray}
we can derive the following variational formulations:
- typically applied in dolfinx tutorials weak form (Approach no. 1)
\int\_\Omega \nabla u \nabla v d \Omega = \int\_\Omega f v d \Omega
- weak formulation with boundary conditions imposed weakly (Approach no. 2)
\begin{eqnarray} \int\_\Omega \nabla u \nabla v d \Omega -\int\_{\Gamma\_D} u \frac{\partial v}{\partial n} d \Gamma\_D -\int\_{\Gamma\_D} \frac{\partial u}{\partial n} v d \Gamma\_D = \int\_\Omega f v d\Omega -\int\_{\Gamma\_D} \bar{u} \frac{\partial v}{\partial n} d \Gamma\_D \end{eqnarray}
- inverse formulation with boundary conditions imposed weakly (Approach no. 3)
\begin{eqnarray} -\int\_\Omega u \nabla^2v d \Omega -\int\_{\Gamma\_D} \frac{\partial u}{\partial n} v d \Gamma\_D = \int\_\Omega f v d\Omega -\int\_{\Gamma\_D} \bar{u} \frac{\partial v}{\partial n} d \Gamma\_D \end{eqnarray}
- original formulation with boundary conditions imposed weakly (Approach no. 4)
\begin{eqnarray} - \int\_\Omega \nabla^2 u v d \Omega - \int\_{\Gamma\_D} u \frac{\partial v}{\partial n} d \Gamma\_D = \int\_\Omega f v d\Omega - \int\_{\Gamma\_D} \bar{u} \frac{\partial v}{\partial n} d \Gamma\_D \end{eqnarray}
The above equations are implemented in the code below. However, while the first two (weak formulations) work seamlessly, inverse and original formulation return NaN solutions.
It seems, that Laplacian cannot be applied on trial or test function, despite the fact that ufl notation allows us to compute Laplacian for the function expressed with the use of the same function space as `u` or `v` (vide line 15: `f = -div(grad(uD))` ).
Is it desired limitation of `DOLFINx`?
Is it possible to force `DOLFINx` to solve such problems?
```auto
from mpi4py import MPI
import dolfinx
from dolfinx.fem.petsc import LinearProblem
from dolfinx import mesh
from ufl import (FacetNormal, SpatialCoordinate, TestFunction, TrialFunction,
div, dot, dx, grad, inner, ds)
import pyvista
domain = mesh.create_unit_square(MPI.COMM_WORLD, 8, 8, mesh.CellType.quadrilateral)
V = dolfinx.fem.functionspace(domain, ("Lagrange", 3))
# manufactured solution: y^2 + 1
u_ex = lambda x: 1+x[1]**2
uD = dolfinx.fem.Function(V)
uD.interpolate( u_ex )
f = -div(grad(uD))
n = FacetNormal(domain)
x = SpatialCoordinate(domain)
u = TrialFunction(V)
v = TestFunction(V)
# ###################################################################
# # # approach #1 - weak standard
# tdim = domain.topology.dim
# fdim = tdim - 1
# domain.topology.create_connectivity(fdim, tdim)
# boundary_facets = mesh.exterior_facet_indices(domain.topology)
# boundary_dofs = fem.locate_dofs_topological(V, fdim, boundary_facets)
# bc = [fem.dirichletbc(uD, boundary_dofs)]
# a = inner(grad(u), grad(v)) * dx
# L = f * v * dx
# # ###################################################################
# # approach #2 - weak w/ BC imposed weakly
bc = []
a = inner(grad(u), grad(v)) * dx \
- inner(u, inner(grad(v), n)) * ds \
- inner(inner(grad(u), n), v) * ds
L = f * v * dx - inner( uD, inner(grad(v), n)) * ds
# ###################################################################
# # # approach #3 - inverse w/ BC imposed weakly
# bc = []
# a = - inner(u, div(grad(v))) * dx \
# - inner(inner(grad(u), n), v) * ds
# L = f * v * dx - inner( uD, inner(grad(v), n)) * ds
# ###################################################################
# # # approach #4 - original w/ BC imposed weakly
# bc = []
# a = - inner(div(grad(u)), v) * dx \
# - inner(u, inner(grad(v), n)) * ds
# L = f * v * dx - inner( uD, inner(grad(v), n)) * ds
problem = LinearProblem(a, L, bcs=bc, petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
uh = problem.solve()
# plotting
pyvista_cells, cell_types, geometry = dolfinx.plot.vtk_mesh(V)
grid = pyvista.UnstructuredGrid(pyvista_cells, cell_types, geometry)
grid.point_data["u"] = uh.x.array
grid.set_active_scalars("u")
plotter = pyvista.Plotter()
plotter.add_text("uh", position="upper_edge", font_size=14, color="black")
plotter.add_mesh(grid, show_edges=True)
plotter.view_xy()
plotter.show()
```
---
### Author:  [@dokken](https://fenicsproject.discourse.group/u/dokken)
#### Post date: [November 12, 2024, 4:22pm UTC](https://fenicsproject.discourse.group/t/is-it-possible-to-solve-original-or-weak-formulations-of-the-boundary-value-problems-in-dolfinx/16282/2 "2024-11-12T16:22:26Z")
> [@marbor](#):
>
> weak formulation with boundary conditions imposed weakly (Approach no. 2)
This is not a stable approach for setting boundary conditions weakly. Consider the excellent introduction in: [https://epubs.siam.org/doi/abs/10.1137/S0036142901384162](https://epubs.siam.org/doi/abs/10.1137/S0036142901384162) (section 2.1).
For the other issues, you need to think about what polynomial spaces you are using for your PDE. Even if you are using a 3rd order Lagrange function space, it is only piecewise continuous (there is a kink between each element. See for instance: [Galerkin methods — FEniCS Workshop](https://jsdokken.com/FEniCS-workshop/src/finite_element_method/galerkin.html#choice-of-basis-functions)
I.e. the gradient is discontinuous, thus, what is a well defined double derivative?
---
### Author:  [@marbor](https://fenicsproject.discourse.group/u/marbor)
#### Post date: [November 12, 2024, 10:16pm UTC](https://fenicsproject.discourse.group/t/is-it-possible-to-solve-original-or-weak-formulations-of-the-boundary-value-problems-in-dolfinx/16282/3 "2024-11-12T22:16:58Z")
Thanks for answering - I’ll give cited paper a look, but I feel like you haven’t answered my question explicitly… 😉 So how is it really: is there a way to solve inverse/original formulations or not?
When it comes to ‘other issues’:
1. I used Lagrange continuous elements just for example purposes, taking care only about differentiability.
Firstly, sometimes continuity of the solution is not necessary. Secondly - if we were interested in continuous solution we could always use Hermitian or spline elements, couldn’t we?
2. Why I’m asking about the issue at all?
I think the possibility of handling original and inverse formulation would make solver more versatile. For example, Trefftz finite element methods are derived from original variational formulation. Similarly, there exist meshless solvers based on original (strong) variational formulation [https://e6.ijs.si/ParallelAndDistributedSystems/publications/69777155.pdf](https://e6.ijs.si/ParallelAndDistributedSystems/publications/69777155.pdf).
3. How can we solve biharmonic equation in DOLFINx if it’s weak form contains \int \nabla^2 u \nabla^2 v~d\Omega?
---
### Author:  [@dokken](https://fenicsproject.discourse.group/u/dokken)
#### Post date: [November 12, 2024, 11:26pm UTC](https://fenicsproject.discourse.group/t/is-it-possible-to-solve-original-or-weak-formulations-of-the-boundary-value-problems-in-dolfinx/16282/4 "2024-11-12T23:26:12Z")
Biharmonic equation: [Biharmonic equation — DOLFINx 0.10.0.0 documentation](https://docs.fenicsproject.org/dolfinx/main/python/demos/demo_biharmonic.html)
Honestly, I can’t give you a full review of Galerkin based finite element methods.
I think a good source of information is: [pde - What is the purpose of using integration by parts in deriving a weak form for FEM discretization? - Computational Science Stack Exchange](https://scicomp.stackexchange.com/questions/7845/what-is-the-purpose-of-using-integration-by-parts-in-deriving-a-weak-form-for-fe#:~:text=bearing%20the%20same%20spirit%20with,not%20require%20integration%20by%20parts).
> [@marbor](#):
>
> - I used Lagrange continuous elements just for example purposes, taking care only about differentiability.
> Firstly, sometimes continuity of the solution is not necessary. Secondly - if we were interested in continuous solution we could always use Hermitian or spline elements, couldn’t we?
Yes, if you want higher order regularity, you would need to use hermitian elements or splines.
> [@marbor](#):
>
> Put me right if I go astray, wouldn’t it be as easy as defining new function spaces/new types of finite element and appropriate operators to generate whole bunch of other numerical methods based on variational formulations if UFL could properly interpret the notation
You say that it is «as easy as». Sure, if someone has time to implement these elements, dofmaps and other constructs, similar to what was done in
> **[GitHub - david-kamensky/tIGAr: A Python library for isogeometric analysis (IGA)...](https://github.com/david-kamensky/tIGAr)**
>
> A Python library for isogeometric analysis (IGA) using FEniCS.
an IGA library built on top of legacy fenics.
However, to make an efficient and maintainable implementation the developer has to been experienced with the fenics code structure, as well as the method they want to implement.
Since we are an open source project mainly driven by people in academic positions, we value any additional suggestions and proposals to improve the framework, and welcome external contributions.
---
### Author:  [@marbor](https://fenicsproject.discourse.group/u/marbor)
#### Post date: [November 21, 2024, 12:36pm UTC](https://fenicsproject.discourse.group/t/is-it-possible-to-solve-original-or-weak-formulations-of-the-boundary-value-problems-in-dolfinx/16282/5 "2024-11-21T12:36:17Z")
I’m terribly sorry for such a long response time – some important matters came up and I couldn’t get back earlier.
Thanks for the link to biharmonic example (I don’t know how it could happen that I couldn’t find it earlier).
Ok. So now I see that we can apply `div grad` operator to a function on the LHS. Putting aside instability of the problems with boundary conditions imposed weakly (which in my opinion shouldn’t ruin the solution, and as a matter of fact, they _do_ give decent results for the weak formulation) can you give any hint why my code does not work for the original and inverse formulations?
When it comes to ‘‘as easy as’’ - well, I should have used quotation marks in my comment… By saying that I didn’t mean to diminish the work put into fenics development. On the contrary, UFL seems to be a great tool that can be a core of a framework for automated generation of many different numerical methods (not necesserily narrowed to FEM-like appproches) derived from variational formulations.
Once again thank you for making me aware of tIGAr library.
And the last question: could you recommend some starting off material for someone interested in understanding how FEniCS/DOLFINx works?
---
### Author:  [@dokken](https://fenicsproject.discourse.group/u/dokken)
#### Post date: [November 21, 2024, 1:25pm UTC](https://fenicsproject.discourse.group/t/is-it-possible-to-solve-original-or-weak-formulations-of-the-boundary-value-problems-in-dolfinx/16282/6 "2024-11-21T13:25:51Z")
I’ve made a quite thorough introduction at: [FEniCS workshop — FEniCS Workshop](https://jsdokken.com/FEniCS-workshop/README.html)
There are various other tutorials out there (that i’ve also made with the help of others, listed at): [Tutorials | Jørgen S. Dokken](https://jsdokken.com/tutorials.html)