I’m trying to change the demo example “Mixed formulation for Poisson equation” given here:
https://fenicsproject.org/docs/dolfin/2019.1.0/python/demos/mixed-poisson/demo_mixed-poisson.py.html
I want to compute flux based on a given pressure gradient. So, I don’t need to solve for u and only need to solve for \sigma.
Suppose I’ve obtained solution for u like this:
from dolfin import *
import matplotlib.pyplot as plt
# Create mesh and define function space
mesh = UnitSquareMesh(32, 32)
V = FunctionSpace(mesh, "Lagrange", 1)
def boundary_bottom(x):
return x[1] < DOLFIN_EPS
def boundary_top(x):
return x[1] > 1.0 - DOLFIN_EPS
# Define boundary condition
bc_top = DirichletBC(V, Constant(0.0), boundary_top)
bc_bottom = DirichletBC(V, Constant(1.0), boundary_bottom)
bcs = [bc_top, bc_bottom]
# Define variational problem
u = TrialFunction(V)
v = TestFunction(V)
f = Expression("0", degree=2)
a = inner(grad(u), grad(v))*dx
L = f*v*dx
# Compute solution
u = Function(V)
solve(a == L, u, bcs)
plot(u)
plt.show()
This gives me simple vertical gradient.
Then, I want to use it to find flux. First, I define mixed space (as far as I know I have to always use mixed space for such problems):
# Define finite elements spaces and build mixed space
BDM = FiniteElement("BDM", mesh.ufl_cell(), 1)
DG = FiniteElement("DG", mesh.ufl_cell(), 0)
W = FunctionSpace(mesh, BDM * DG)
# Define trial and test functions
(sigma, u) = TrialFunctions(W)
(tau, v) = TestFunctions(W)
And now I want to substitute u for the u I’ve obtained later. And use it later to find the flux in:
\begin{split}\sigma - \nabla u &= 0 \quad {\rm in} \ \Omega, \\ \nabla \cdot \sigma &= - f \quad {\rm in} \ \Omega,\end{split}
I’ve tried different way of changing the demo example, but got different errors.
Something tells me that my approach is wrong from the very beginning, from the step of defining mixed space or test and trial function.