Hello everyone concerned,
Recently, I read one paper whose DOI is 10.1016/j.cam.2014.04.022.
The algorithm provided is eq.(10) on page 425.
I try to implement this algorithm by FEniCS.
However, I wonder how to implement the third term on RHS.
This term is not a inner product on the Dirichlet boundary.
It’s just a scalar product I think.
Moreover, as the functional space defines, the test function on this boundary is a constant without an exact value.
Does someone knows that how to implement this term by FEniCS?
Thanks very much.
Here is my code:
from fenics import *
import matplotlib.pyplot as plt
import numpy as np
import math
def g(u):
return u**2+u
T = 1.0
num_steps = 2**3
dt = T / num_steps
t = 0.0
nx = 2**3
# Define exact solution
u_e = Expression(('t*(cos(x[0])+x[0])+x[0]'), t = 0.0, degree = 1)
# u_e = Expression("1 + x[0]*x[0] + 2*x[1]*x[1]", degree = 2)
# Create mesh and define function space
mesh = IntervalMesh(nx, 0, pi/2)
V = FunctionSpace(mesh, 'P', 1)
plot(mesh)
plt.savefig('mesh.png')
# Define boundary condition
class Left(SubDomain):
def inside(self, x, on_boundary):
return near(x[0], 0.0)
class Right(SubDomain):
def inside(self, x, on_boundary):
return near(x[0], pi/2)
# Initialize sub-domain instances
left = Left()
right = Right()
# Initialize mesh function for interior domains
domains = MeshFunction("size_t", mesh, mesh.topology().dim())
domains.set_all(0)
# Initialize mesh function for boundary domains
boundaries = MeshFunction("size_t", mesh, mesh.topology().dim()-1)
boundaries.set_all(0)
left.mark(boundaries, 1)
right.mark(boundaries, 2)
# Define input data
E = Expression(('-2*t-t+1/24*pow(pi,3)+5/4*pow(t,2)*pi+1/8*t*pi*pi\
+t*pi+1/12*t*pow(pi,3)+1/24*pow(t,2)*pow(pi,3)+1/8*pow(pi,2)'),\
t = 0.0, degree = 0)
dE = Expression(('-2-1+5/2*t*pi+1/8*pi*pi+pi+1/12*pi*pi*pi \
+1/12*t*pi*pi*pi'), t = 0.0, degree = 0)
# f = Expression(('-pow((cos(x[0])+x[0]+1),2)*pow(t,3)/3' \
# + '-(cos(x[0])+x[0]+1)*x[0]*pow(t,2)' \
# + '+(2*pow((cos(x[0])+x[0]),2) + cos(x[0]) - pow(x[0],2))*t' \
# + '+(2*x[0] + 1)*(cos(x[0])+x[0])'), t = 0.0, degree = 2)
f = Expression(('(2*pow((cos(x[0])+x[0]),2) + cos(x[0]) )*t' \
+ '+(2*x[0] + 1)*(cos(x[0])+x[0])'), t = 0.0, degree = 2)
h = Expression(('t + 1'), t = 0.0, degree=0)
# exterior boundaries
dx = Measure('dx', domain=mesh, subdomain_data=domains)
ds = Measure('ds', domain=mesh, subdomain_data=boundaries)
u_1 = interpolate(u_e,V)
du = TrialFunction(V)
v = TestFunction(V)
u = Function(V)
F = inner(pow(u,2),v)*dx(0) + inner(u,v)*dx(0) + \
dt*inner(grad(u), grad(v))*dx(0) -\
dt*inner(f,v)*dx(0) + dt*inner(h,v)*ds(1) \
-dt*dot(v,dE - inner(f,1)*dx(0) + inner(h,1)*ds(1))*ds(2) \
-inner(pow(u_1,2),v)*dx(0) - inner(u_1,v)*dx(0)
J = derivative(F, u, du)
problem = NonlinearVariationalProblem(F, u, J)
solver = NonlinearVariationalSolver(problem)
t = dt
# Create VTK file for saving solution
file1 = File('solution/u/u.pvd')
file2 = File('solution/ue/ue.pvd')
# while t<T:
# u_e.t = t
# bcs.t = t
# f.t = t
# u0 = interpolate(u_e,V)
# problem = NonlinearVariationalProblem(F, u, bcs, J)
# solver = NonlinearVariationalSolver(problem)
# solver.solve()
# file1 << (u, t)
# file2 << (u0, t)
# t += dt
# u_1.assign(u)
# maxdiff = np.linalg.norm(np.array(u0.vector())-np.array(u.vector()), ord=2)/np.linalg.norm(np.array(u0.vector()), ord = 2)