#======================================================
#  Forcing function and exact solutions for 1D problem
#======================================================
cA_expression = Expression( \
'exp(x[0])*sin(pi*x[0])*pow(sin(t),2)' \
,t=t)

cB_expression = Expression( \
'1 - sin(5*pi*x[0])*pow(cos(3*pi*t*x[0]),2)*(x[0] - 1) - x[0]' \
,t=t)

cC_expression = Expression( \
'-k*exp(x[0])*sin(pi*x[0])*pow(sin(t),2)*(x[0] + sin(5*pi*x[0])*pow(cos(3*pi*t*x[0]),2)*(x[0] - 1) - 1)' \
,t=t,k=keq)

psiA1_expression = Expression( \
'exp(x[0])*sin(pi*x[0])*pow(sin(t),2) - k*exp(x[0])*sin(pi*x[0])*pow(sin(t),2)*(x[0] + sin(5*pi*x[0])*pow(cos(3*pi*t*x[0]),2)*(x[0] - 1) - 1)' \
,t=t,k=keq)

psiB1_expression = Expression( \
'1 - sin(5*pi*x[0])*pow(cos(3*pi*t*x[0]),2)*(x[0] - 1) - k*exp(x[0])*sin(pi*x[0])*pow(sin(t),2)*(x[0] + sin(5*pi*x[0])*pow(cos(3*pi*t*x[0]),2)*(x[0] - 1) - 1) - x[0]' \
,t=t,k=keq)

psiA2_expression = Expression( \
'x[0] + exp(x[0])*sin(pi*x[0])*pow(sin(t),2) + sin(5*pi*x[0])*pow(cos(3*pi*t*x[0]),2)*(x[0] - 1) - 1' \
,t=t)

psiC2_expression = Expression( \
'1 - sin(5*pi*x[0])*pow(cos(3*pi*t*x[0]),2)*(x[0] - 1) - k*exp(x[0])*sin(pi*x[0])*pow(sin(t),2)*(x[0] + sin(5*pi*x[0])*pow(cos(3*pi*t*x[0]),2)*(x[0] - 1) - 1) - x[0]' \
,t=t,k=keq)

fpsiA1_expression = Expression( \
'D*(pow(pi,2)*exp(x[0])*sin(pi*x[0])*pow(sin(t),2) - 2*pi*exp(x[0])*cos(pi*x[0])*pow(sin(t),2) - exp(x[0])*sin(pi*x[0])*pow(sin(t),2) + 2*k*exp(x[0])*sin(pi*x[0])*pow(sin(t),2)*(sin(5*pi*x[0])*pow(cos(3*pi*t*x[0]),2) + 5*pi*cos(5*pi*x[0])*pow(cos(3*pi*t*x[0]),2)*(x[0] - 1) - 6*t*pi*sin(3*pi*t*x[0])*sin(5*pi*x[0])*cos(3*pi*t*x[0])*(x[0] - 1) + 1) + k*exp(x[0])*sin(pi*x[0])*pow(sin(t),2)*(x[0] + sin(5*pi*x[0])*pow(cos(3*pi*t*x[0]),2)*(x[0] - 1) - 1) - k*exp(x[0])*sin(pi*x[0])*pow(sin(t),2)*(25*pow(pi,2)*sin(5*pi*x[0])*pow(cos(3*pi*t*x[0]),2)*(x[0] - 1) - 10*pi*cos(5*pi*x[0])*pow(cos(3*pi*t*x[0]),2) + 18*pow(t,2)*pow(pi,2)*sin(5*pi*x[0])*pow(cos(3*pi*t*x[0]),2)*(x[0] - 1) - 18*pow(t,2)*pow(pi,2)*pow(sin(3*pi*t*x[0]),2)*sin(5*pi*x[0])*(x[0] - 1) + 12*t*pi*sin(3*pi*t*x[0])*sin(5*pi*x[0])*cos(3*pi*t*x[0]) + 60*t*pow(pi,2)*sin(3*pi*t*x[0])*cos(5*pi*x[0])*cos(3*pi*t*x[0])*(x[0] - 1)) - k*pow(pi,2)*exp(x[0])*sin(pi*x[0])*pow(sin(t),2)*(x[0] + sin(5*pi*x[0])*pow(cos(3*pi*t*x[0]),2)*(x[0] - 1) - 1) + 2*k*pi*exp(x[0])*cos(pi*x[0])*pow(sin(t),2)*(sin(5*pi*x[0])*pow(cos(3*pi*t*x[0]),2) + 5*pi*cos(5*pi*x[0])*pow(cos(3*pi*t*x[0]),2)*(x[0] - 1) - 6*t*pi*sin(3*pi*t*x[0])*sin(5*pi*x[0])*cos(3*pi*t*x[0])*(x[0] - 1) + 1) + 2*k*pi*exp(x[0])*cos(pi*x[0])*pow(sin(t),2)*(x[0] + sin(5*pi*x[0])*pow(cos(3*pi*t*x[0]),2)*(x[0] - 1) - 1)) - v*(k*exp(x[0])*sin(pi*x[0])*pow(sin(t),2)*(sin(5*pi*x[0])*pow(cos(3*pi*t*x[0]),2) + 5*pi*cos(5*pi*x[0])*pow(cos(3*pi*t*x[0]),2)*(x[0] - 1) - 6*t*pi*sin(3*pi*t*x[0])*sin(5*pi*x[0])*cos(3*pi*t*x[0])*(x[0] - 1) + 1) - pi*exp(x[0])*cos(pi*x[0])*pow(sin(t),2) - exp(x[0])*sin(pi*x[0])*pow(sin(t),2) + k*exp(x[0])*sin(pi*x[0])*pow(sin(t),2)*(x[0] + sin(5*pi*x[0])*pow(cos(3*pi*t*x[0]),2)*(x[0] - 1) - 1) + k*pi*exp(x[0])*cos(pi*x[0])*pow(sin(t),2)*(x[0] + sin(5*pi*x[0])*pow(cos(3*pi*t*x[0]),2)*(x[0] - 1) - 1)) + 2*exp(x[0])*sin(pi*x[0])*cos(t)*sin(t) - 2*k*exp(x[0])*sin(pi*x[0])*cos(t)*sin(t)*(x[0] + sin(5*pi*x[0])*pow(cos(3*pi*t*x[0]),2)*(x[0] - 1) - 1) + 6*k*x[0]*pi*sin(3*pi*t*x[0])*exp(x[0])*sin(pi*x[0])*sin(5*pi*x[0])*pow(sin(t),2)*cos(3*pi*t*x[0])*(x[0] - 1)' \
,t=t,D=diffusivity,v=vel,k=keq)

fpsiB1_expression = Expression( \
'D*(10*pi*cos(5*pi*x[0])*pow(cos(3*pi*t*x[0]),2) - 25*pow(pi,2)*sin(5*pi*x[0])*pow(cos(3*pi*t*x[0]),2)*(x[0] - 1) + 2*k*exp(x[0])*sin(pi*x[0])*pow(sin(t),2)*(sin(5*pi*x[0])*pow(cos(3*pi*t*x[0]),2) + 5*pi*cos(5*pi*x[0])*pow(cos(3*pi*t*x[0]),2)*(x[0] - 1) - 6*t*pi*sin(3*pi*t*x[0])*sin(5*pi*x[0])*cos(3*pi*t*x[0])*(x[0] - 1) + 1) + k*exp(x[0])*sin(pi*x[0])*pow(sin(t),2)*(x[0] + sin(5*pi*x[0])*pow(cos(3*pi*t*x[0]),2)*(x[0] - 1) - 1) - k*exp(x[0])*sin(pi*x[0])*pow(sin(t),2)*(25*pow(pi,2)*sin(5*pi*x[0])*pow(cos(3*pi*t*x[0]),2)*(x[0] - 1) - 10*pi*cos(5*pi*x[0])*pow(cos(3*pi*t*x[0]),2) + 18*pow(t,2)*pow(pi,2)*sin(5*pi*x[0])*pow(cos(3*pi*t*x[0]),2)*(x[0] - 1) - 18*pow(t,2)*pow(pi,2)*pow(sin(3*pi*t*x[0]),2)*sin(5*pi*x[0])*(x[0] - 1) + 12*t*pi*sin(3*pi*t*x[0])*sin(5*pi*x[0])*cos(3*pi*t*x[0]) + 60*t*pow(pi,2)*sin(3*pi*t*x[0])*cos(5*pi*x[0])*cos(3*pi*t*x[0])*(x[0] - 1)) - 18*pow(t,2)*pow(pi,2)*sin(5*pi*x[0])*pow(cos(3*pi*t*x[0]),2)*(x[0] - 1) + 18*pow(t,2)*pow(pi,2)*pow(sin(3*pi*t*x[0]),2)*sin(5*pi*x[0])*(x[0] - 1) - 12*t*pi*sin(3*pi*t*x[0])*sin(5*pi*x[0])*cos(3*pi*t*x[0]) - k*pow(pi,2)*exp(x[0])*sin(pi*x[0])*pow(sin(t),2)*(x[0] + sin(5*pi*x[0])*pow(cos(3*pi*t*x[0]),2)*(x[0] - 1) - 1) + 2*k*pi*exp(x[0])*cos(pi*x[0])*pow(sin(t),2)*(sin(5*pi*x[0])*pow(cos(3*pi*t*x[0]),2) + 5*pi*cos(5*pi*x[0])*pow(cos(3*pi*t*x[0]),2)*(x[0] - 1) - 6*t*pi*sin(3*pi*t*x[0])*sin(5*pi*x[0])*cos(3*pi*t*x[0])*(x[0] - 1) + 1) + 2*k*pi*exp(x[0])*cos(pi*x[0])*pow(sin(t),2)*(x[0] + sin(5*pi*x[0])*pow(cos(3*pi*t*x[0]),2)*(x[0] - 1) - 1) - 60*t*pow(pi,2)*sin(3*pi*t*x[0])*cos(5*pi*x[0])*cos(3*pi*t*x[0])*(x[0] - 1)) - v*(sin(5*pi*x[0])*pow(cos(3*pi*t*x[0]),2) + 5*pi*cos(5*pi*x[0])*pow(cos(3*pi*t*x[0]),2)*(x[0] - 1) + k*exp(x[0])*sin(pi*x[0])*pow(sin(t),2)*(sin(5*pi*x[0])*pow(cos(3*pi*t*x[0]),2) + 5*pi*cos(5*pi*x[0])*pow(cos(3*pi*t*x[0]),2)*(x[0] - 1) - 6*t*pi*sin(3*pi*t*x[0])*sin(5*pi*x[0])*cos(3*pi*t*x[0])*(x[0] - 1) + 1) + k*exp(x[0])*sin(pi*x[0])*pow(sin(t),2)*(x[0] + sin(5*pi*x[0])*pow(cos(3*pi*t*x[0]),2)*(x[0] - 1) - 1) - 6*t*pi*sin(3*pi*t*x[0])*sin(5*pi*x[0])*cos(3*pi*t*x[0])*(x[0] - 1) + k*pi*exp(x[0])*cos(pi*x[0])*pow(sin(t),2)*(x[0] + sin(5*pi*x[0])*pow(cos(3*pi*t*x[0]),2)*(x[0] - 1) - 1) + 1) - 2*k*exp(x[0])*sin(pi*x[0])*cos(t)*sin(t)*(x[0] + sin(5*pi*x[0])*pow(cos(3*pi*t*x[0]),2)*(x[0] - 1) - 1) + 6*x[0]*pi*sin(3*pi*t*x[0])*sin(5*pi*x[0])*cos(3*pi*t*x[0])*(x[0] - 1) + 6*k*x[0]*pi*sin(3*pi*t*x[0])*exp(x[0])*sin(pi*x[0])*sin(5*pi*x[0])*pow(sin(t),2)*cos(3*pi*t*x[0])*(x[0] - 1)' \
,t=t,D=diffusivity,v=vel,k=keq)

fpsiA2_expression = Expression( \
'D*(pow(pi,2)*exp(x[0])*sin(pi*x[0])*pow(sin(t),2) - 10*pi*cos(5*pi*x[0])*pow(cos(3*pi*t*x[0]),2) - 2*pi*exp(x[0])*cos(pi*x[0])*pow(sin(t),2) - exp(x[0])*sin(pi*x[0])*pow(sin(t),2) + 25*pow(pi,2)*sin(5*pi*x[0])*pow(cos(3*pi*t*x[0]),2)*(x[0] - 1) + 18*pow(t,2)*pow(pi,2)*sin(5*pi*x[0])*pow(cos(3*pi*t*x[0]),2)*(x[0] - 1) - 18*pow(t,2)*pow(pi,2)*pow(sin(3*pi*t*x[0]),2)*sin(5*pi*x[0])*(x[0] - 1) + 12*t*pi*sin(3*pi*t*x[0])*sin(5*pi*x[0])*cos(3*pi*t*x[0]) + 60*t*pow(pi,2)*sin(3*pi*t*x[0])*cos(5*pi*x[0])*cos(3*pi*t*x[0])*(x[0] - 1)) + v*(sin(5*pi*x[0])*pow(cos(3*pi*t*x[0]),2) + exp(x[0])*sin(pi*x[0])*pow(sin(t),2) + pi*exp(x[0])*cos(pi*x[0])*pow(sin(t),2) + 5*pi*cos(5*pi*x[0])*pow(cos(3*pi*t*x[0]),2)*(x[0] - 1) - 6*t*pi*sin(3*pi*t*x[0])*sin(5*pi*x[0])*cos(3*pi*t*x[0])*(x[0] - 1) + 1) + 2*exp(x[0])*sin(pi*x[0])*cos(t)*sin(t) - 6*x[0]*pi*sin(3*pi*t*x[0])*sin(5*pi*x[0])*cos(3*pi*t*x[0])*(x[0] - 1)' \
,t=t,D=diffusivity,v=vel,k=keq)