Thank you dokken, another boundaries question, I learn another way to set the mixed boundary conditions(22. Poisson equation with multiple subdomains — FEniCS Project), namely

for x in [0,Length]

u=(Length, t)=0

u_x(0, t)=func(t)

The first (previous) version is

```
from dolfin import *
import numpy as np
E=3300000000.0
rho=1300.0
c = sqrt(E/rho)
Po=30000000.0
t1=0.00015
t2=0.001
b0=1.0
Length=30.0
# Time variables
dt = 0.0001; t = 0; T = 0.1
mesh = IntervalMesh(1000,0,Length)
V=FunctionSpace(mesh, "Lagrange", 1)
class Right(SubDomain):
def inside(self,x,on_boundary):
return on_boundary and near(x[0], Length)
boundaries = MeshFunction("size_t", mesh,0)
boundaries.set_all(0)
# Now mark your Neumann boundary
rightbound = Right()
rightbound.mark(boundaries, 1)
ds=Measure('ds',subdomain_data=boundaries)
# Previous and current solution
u1= interpolate(Constant(0.0), V)
u0= interpolate(Constant(0.0), V)
# Variational problem at each time
u = TrialFunction(V)
v = TestFunction(V)
blast=Constant((Po*(1-t/t2)*np.exp(-b0*((t-t1)/t2)))/E)
def blast_func(t):
return (Po*(1-t/t2)*np.exp(-b0*((t-t1)/t2)))/E
a = u*v*dx + dt*dt*c*c*inner(grad(u), grad(v))*dx
L = 2*u1*v*dx - u0*v*dx + blast*v*ds
bc = DirichletBC(V, 0, rightbound)
u=Function(V)
vtkfile1=File('wave/u.pvd')
vtkfile2=File('wave/v.pvd')
while t <= T:
A, b = assemble_system(a, L, bc)
solve(A, u.vector(), b)
velocity = (u - u0) / dt
u0.assign(u1)
u1.assign(u)
t += dt
blast.assign(blast_func(t))
vtkfile1 << (u,t)
vtkfile2 << (project(velocity, V), t)
```

The second version is

```
from dolfin import *
import numpy as np
E=3300000000.0
rho=1300.0
c = sqrt(E/rho)
Po=30000000.0
t1=0.00015
t2=0.001
b0=1.0
Length=30.0
# Time variables
dt = 0.0001; t = 0; T = 0.1
mesh = IntervalMesh(1000,0,Length)
V=FunctionSpace(mesh, "Lagrange", 1)
class LeftBoundary(SubDomain):
def inside(self, x, on_boundary):
return near (x[0], 0.0)
class RightBoundary(SubDomain):
def inside(self, x, on_boundary):
return near (x[0], Length)
RB=RightBoundary()
LB=LeftBoundary()
boundaries = MeshFunction("size_t", mesh, 0)
boundaries.set_all(0)
RB.mark(boundaries, 1)
LB.mark(boundaries, 2)
blast=Constant((Po*(1-t/t2)*np.exp(-b0*((t-t1)/t2)))/E)
def blast_func(t):
return (Po*(1-t/t2)*np.exp(-b0*((t-t1)/t2)))/E
# Previous and current solution
u1= interpolate(Constant(0.0), V)
u0= interpolate(Constant(0.0), V)
# Variational problem at each time
u = TrialFunction(V)
v = TestFunction(V)
a = u*v*dx + dt*dt*c*c*inner(grad(u), grad(v))*dx
L = 2*u1*v*dx - u0*v*dx + blast*v*ds(2)
bc = DirichletBC(V, 0.0, boundaries, 1)
u=Function(V)
vtkfile1=File('wave/u.pvd')
vtkfile2=File('wave/v.pvd')
while t <= T:
A, b = assemble_system(a, L, bc)
solve(A, u.vector(), b)
velocity = (u - u0) / dt
u0.assign(u1)
u1.assign(u)
t += dt
blast.assign(blast_func(t))
vtkfile1 << (u,t)
vtkfile2 << (project(velocity, V), t)
```

In the first one, the wave will rebound from the right side, the second is all zero, I don’t know where the problem is in my second version code.