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?
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()
So how is it really: is there a way to solve inverse/original formulations or not?