How to apply boundary conditions and solve different trial and test function space system?

Question

, where Omega is a 3D UnitCubeMesh.

I tried the below code, but I got a singular system Au=b and obtained a wrong solution.

from fenics import *
import numpy as np 
from numpy import linalg

mesh = UnitCubeMesh(10, 10, 10)
RT = FunctionSpace(mesh, 'RT', 1)
DG = FunctionSpace(mesh, 'DG', 0)
Nedelec = FunctionSpace(mesh, 'N1curl', 1)

u = TrialFunction(RT)
q = TestFunction(DG)
eta = TestFunction(Nedelec)

bc_RT = DirichletBC(RT, Constant([0.0, 0.0, 0.0]), DomainBoundary())
bc_Nedelec = DirichletBC(Nedelec, Constant([0.0, 0.0, 0.0]), DomainBoundary())

a = div(u)*q*dx
L = phi*q*dx
a_ = inner(u, curl(eta))*dx
dummy = inner(Constant([0,0,0]), eta)*dx

A, b = assemble_system(a, L, bc_RT)
A_, _ = assemble_system(a_, dummy, [bc_RT, bc_Nedelec])

AA = np.vstack((A.array(), A_.array())) # merge two blocks of matrix(stiffness matrix)
zero = np.zeros(A_.array().shape[0])
bb = np.hstack((b.array(), zero)) # merge two blocks of vector(source term)

uu = linalg.lstsq(AA, bb)[0] # use least square algorithm in numpy

u = Function(RT)
u.vector().set_local(uu)

plot(u, title='solution u')

Comments

Hi, what is the strong form of this problem?

---

Thanks for your answer!!!

I didn't get a strong form from my supervisor. He told me the above linear form is solvable and only exists one solution.
What's more, the Raviart-Thomas element space for H0(div, Omega) of lowest order; the piecewise constant element space for L_0_2(Omega); the continuous linear element space for H_0_1(Omega).

Sorry about the description, my native language is not English.

---

Okay, the reason I asked is that yours resembles the problem of finding a vector $u$ field given its divergence $\phi$. This is the first eq. in your weak form. Observing that $\nabla\cdot\nabla\times$ is zero you will need to constrain $u$ to get a well posed problem. That constraint is the second eq. in your weak form. Alternatively you could make an ansatz $u=\nabla\psi$ for some $\psi$ to be found. This way the curl constraint is automatically satisfied. See also here.

---

Thanks very much!

Actually, I understand your explanations.
Before applying boundary condition, I get Au=b from eq1, A_u=0. Both A and A_ are not square linear systems.
Although from theory the curl constraint is automatically satisfied, I also need explicitly apply boundary condition to A and A_, right?
And then I merge A and A_ to form AA, merge b and 0 to form bb.
Finally I use least square method to solve AA*u=bb.

The above process is right?

---

With the ansatz I suggested you would end up with a Laplace problem for $\psi$. About applying boundary conditions; if the matrix is square you want to zero row $i$ corresponding to some boundary degree of freedom and put 1 on the diagonal, agree? For a rectangular matrix if $V$ is the space of test functions and you have DirichletBC(V, ...) zeroing rows is not a problem. But what is the step equivalent to setting the diagonal?

---

Thank you for your patience!!!

I update the source code.

My supervisor told me not replacing the original problem with a Laplace problem for ψ.

I agree "if the matrix is square you want to zero row i corresponding to some boundary degree of freedom and put 1 on the diagonal".

I notice the functions zero() and zero_columns() here.
And I do some numerical experiments. I find that when the linear system is not symmetric, applying DirichletBC() means: zeroing row i corresponding to some boundary DOFs for test functions, zeroing column j corresponding to some boundary DOFs for trial functions(at the same time updating a (right-hand side) vector to reflect the changes).

---

Yes, using zero and zero_columns might be a way to go if you must stick with the rectangular system.