This is a read only copy of the old FEniCS QA forum. Please visit the new QA forum to ask questions

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

+1 vote

enter image description here

, 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')
asked Jul 7, 2017 by fanronghong FEniCS User (1,680 points)
edited Jul 9, 2017 by fanronghong

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.

...