Sure. Here is the complete code. By global LHS I mean the a matrix which is constructed from a_00, a_01 and a_10 in the code. the none is supposed to be zero matrix.
and by global RHS I mean the L vector.
from __future__ import print_function
from fenics import *
import dolfin
import numpy as np
import matplotlib
from matplotlib import pyplot
import matplotlib.pyplot as plt
T = 0.1 # final time
num_steps = 50 # number of time steps
dt = T / num_steps # time step size
mu = 1.0 # kinematic viscosity
rho = 1 # density
# Create mesh and define function spaces
P = 100
mesh = UnitSquareMesh(P, P)
u_0 = Expression(('-cos(2*pi*x[0])*sin(2*pi*x[1])', 'cos(2*pi*x[1])*sin(2*pi*x[0])'), degree = 1)
#u_0 = Expression(('sin(2*pi*x[1])', 'cos(2*pi*x[0])'), degree = 1)
class PeriodicBoundary(SubDomain):
# Left boundary is "target domain" G
def inside(self, x, on_boundary):
return (near(x[0], 0) or near(x[1],0)) and on_boundary
# Map right boundary (H) to left boundary (G)
def map(self, x, y):
if near(x[0], 1):
y[0] = x[0] - 1
y[1] = x[1]
elif near(x[1], 1):
y[0] = x[0]
y[1] = x[1] - 1
else:
y[0] = x[0]
y[1] = x[1]
# Define boundary condition
pbc = PeriodicBoundary()
V = FunctionSpace(mesh, 'RT', 2, constrained_domain=pbc)
Q = FunctionSpace(mesh, 'DG', 1,constrained_domain=pbc)
# # Define boundaries
# inflow = 'near(x[0], 0)'
# outflow = 'near(x[0], 1)'
# walls = 'near(x[1], 0) || near(x[1], 1)'
pressure_PC = '(near(x[0],0) and near(x[1],0)) || (near(x[0],0) and near(x[1],1)) || (near(x[0],1) and near(x[1],0)) || (near(x[0],1) and near(x[1],1))'
# # Define boundary conditions
# bcu_noslip = DirichletBC(V, Constant((0, 0)), walls)
# bcp_inflow = DirichletBC(Q, Constant(8), inflow)
# bcp_outflow = DirichletBC(Q, Constant(0), outflow)
bcp = DirichletBC(Q, Constant(0), pressure_PC)
# bcu = [bcu_noslip]
bcp = [bcp]
# Define trial and test functions
u = TrialFunction(V)
v = TestFunction(V)
p = TrialFunction(Q)
q = TestFunction(Q)
# Define functions for solutions at previous and current time steps
u_n = Function(V)
u_nn = Function(V)
u_nn = project(u_0, V)
u_n = project(u_0, V)
u_ = Function(V)
p_n = Function(Q)
p_ = Function(Q)
# Define expressions used in variational forms
n = FacetNormal(mesh)
f = Constant((0, 0))
k = Constant(dt)
mu = Constant(mu)
rho = Constant(rho)
# Define variational problem for step 2
####################
a_01 = -inner(p,div(v))*dx
a_10 = -inner(div(u),q)*dx
L_0 = inner(f,v)*dx
L_1 = inner(Constant(1.0), q)*dx
def jump(phi, n):
return outer(phi("+"), n("+")) + outer(phi("-"), n("-"))
# Time-stepping
t = 0
for n in range(num_steps):
# Update current time
t += dt
lmbda = conditional(gt(dot(u_n, n), 0), 1, 0)
u_uw = lmbda("+") * u("+") + lmbda("-") * u("-")
a_00 =(inner(u / k, v) * dx
- inner(u, div(outer(v, u_n))) * dx
+ inner((dot(u_n, n))("+") * u_uw, v("+")) * dS
+ inner((dot(u_n, n))("-") * u_uw, v("-")) * dS
+ inner(dot(u_n, n) * lmbda * u, v) * ds)
a_00 = assemble(a_00)
a_01 = assemble(a_01)
a_10 = assemble(a_10)
none_matrix = np.zeros_like(a_00) * np.nan
a = [[a_00,a_01], [a_10,none_matrix]]
a = assemble(a)
L_0 += inner(u_n/k , v) * dx - inner(dot(u_n,n)*(1-lmbda)*u_0,v) * ds
L_0 = assemble(L_0)
L_1 = assemble(L_1)
L = ([L_0,L_1])
L = assemble(L)
solve(a, u_.vector(), L)
# Plot solution
plot(u_)
print('max u:', u_.vector().max())