Different boundary conditions for trial and test function space

Question

I would like to know if there is an easy way to implement a method with different boundary conditions on the trial and test function spaces. Thanks.

The simplest example is a second-order ODE:
$$\ddot u = f, \qquad u(0)=u_0, \quad \dot u(0)=u_1.$$
A simple weak formulation at the continuous level is: find $u\in H^1([0,T])$ with $u(0)=u_0$ such that
$$ -\int_0^T \dot u \dot v = u_1v(0)+\int_0^T fv=\int_0^T (f+u_1\delta)v+, \qquad \forall v\in \{v\in H^1([0,T])|v(T)=0\}.$$
The above can be discretized using, for example CG1 for both u and v. The implementation difficulty is that the boundary condition for the trial function u is $u(0)=u_0$ but the boundary condition for the test function v is $v(T)=0$. In effect, Dirichlet boundary condition for u is applied as usual, but the equation corresponding to testing against the basis supported at t=T is dropped.

Update:
Below is a "dirty implementation" using FEniCS and numpy:

from dolfin import *
import numpy as np
# input data
f = Expression("-225.0*sin(15.0*x[0])")
u0 = 0.0
u1 = 15.0
u_ext = Expression("sin(15.0*x[0])")
# discretization parameters
mesh = UnitIntervalMesh(300)
V = FunctionSpace(mesh, 'CG', 1)
# galerkin method
u = TrialFunction(V)
v = TestFunction(V)
a = -(u.dx(0)*v.dx(0))*dx
L = f*v*dx
# assemble and apply PointSource on rhs
A = assemble(a)
b = assemble(L)
PointSource(V, Point(0.0), u1).apply(b)

#--- dirty part begins ---#
## export to numpy
A = A.array()
b = b.array()
## remove the equation corresponding to testing
#  against the test function centered at t=1.0.
## for some reason, FEniCS CG1 dofs on the unit
#  interval is numbered in reverse. See,
#      V.dofmap.tabulate_all_coordinates(mesh)
#  Thus we want to remove the very first equation
A = A[1:,] 
b = b[1:]
# add the equation for the dirichlet BC: u(0)=u_0
row = np.zeros(A.shape[1])
row[-1] = 1.0 # the last component corresponds to t=0.0.
A = np.vstack([A, row])
b = np.hstack([b, u0])
# numpy solve.
x = np.linalg.solve(A, b)
# plug the solution back to FEniCS
uh = Function(V)
uh.vector()[:] = x
#--- dirty part ends ---#

# compute the L2 error
Vb = FunctionSpace(mesh, 'CG', 2)
e = interpolate(u_ext, Vb) - interpolate(uh, Vb)
print("L2 error={}".format(np.sqrt(assemble(e*e*dx))))
#plot(uh, interactive=True)

The method itself is basically the standard central difference written in a finite element way (with better handling of the right-hand side). Similar time stepping can be used for the linear wave equation as well. But the dirty implementation here is far too dirty for PDE problems.

Comments

Is your question whether or not the 'Dirichlet' dofs can be eliminated from the matrix?

---

Update: More clarification and a FEniCS+numpy implementation is added.
It can probably be interpreted as removing Dirchlet dof from the matrix. More simply, it is to shrink the test function space.

---

I am wondering if Restriction can be used to do this?

Answer 1

Larry, You can do this by not imposing any Dirichlet boundary condition directly but adding an addition real number parameter in both the trial and the test space. So the variational formulation becomes: find $u\in H^1$, $r\in \mathbb R$ such that

$$-\int \dot u\dot v + u(0) s + r v(T) = \int fv, \quad v\in H^1,\
s\in \mathbb R.$$

(I am assuming $u_0=u_1=0$ for simplicity only.) Here is an implementation:

from dolfin import *
import numpy as np
# problem is u'' = f, u(0) = u'(0) = 0  with  u = cos(15 x) -1
f = Expression("-225.0*cos(15.0*x[0])")
mesh = UnitIntervalMesh(300)
V = FunctionSpace(mesh, 'CG', 1)
R = FunctionSpace(mesh, 'Real', 0)
W = MixedFunctionSpace([V, R])
left = Expression('1-x[0]') # 1 at left endpoint, 0 at right
right = Expression('x[0]')
# Galerkin method
u, r = TrialFunctions(W)
v, s = TestFunctions(W)
a = -(u.dx(0)*v.dx(0))*dx + left*u*s*ds + right*r*v*ds
L = f*v*dx
w = Function(W)
solve(a == L, w)
u, r = w.split()