# -*- coding: utf-8 -*-
"""
Created on Tue Oct  7 14:28:44 2014

@author: hans-werner
"""

from dolfin import *

# -----------------------------------------------------------------------------
# Problem a):
# -----------------------------------------------------------------------------
# Mesh and Function Space
mesh = IntervalMesh(20,-2,1)
V    = FunctionSpace(mesh,"CG",1)

# Define Parameters
fa = Expression('x[0]')

# Define Dirichlet boundary conditions at the left endpoint
def q2a_boundary_left(x,on_boundary):
  return on_boundary and abs(x+2) < DOLFIN_EPS

bc = DirichletBC(V,Constant('2.0'),q2a_boundary_left)

# Formulate the weak form a(u,v) = L(v)
u = TrialFunction(V)
v = TestFunction(V)
a = inner(grad(u),grad(v))*dx + u*v*dx
L = fa*v*dx + Constant('5.0')*v*ds

# Solve the problem
u = Function(V)
solve(a==L,u,bc)

# Make a plot
plot(u,interactive=True)


# -----------------------------------------------------------------------------
# Problem b)
# -----------------------------------------------------------------------------
# Mesh and Function Space
mesh = IntervalMesh(20,0.0,1.0)
V    = FunctionSpace(mesh,"CG",1)

# Forcing term is the same as for a)

# Define the Dirichlet boundary condition on the left
def q2b_boundary_left(x,on_boundary):
  return on_boundary and abs(x) < DOLFIN_EPS
  
bc = DirichletBC(V,Constant('0.0'),q2b_boundary_left)

# Formulate the weak form a(u,v) = L(v)
u = TrialFunction(V)
v = TestFunction(V)
a = inner(grad(u),grad(v))*dx + u*v*dx + 2.0*u.dx(0)*v*dx
L = fa*v*dx

# Solve the problem
u = Function(V)
solve(a==L,u,bc)

# Make a plot
plot(u,interactive=True)


# -----------------------------------------------------------------------------
# Problem c) 
# -----------------------------------------------------------------------------
# Mesh and Function Space is the same as before

# Define diffusion coefficient
q = Expression("1+2*x[0]")

# Define Dirichlet Boundary Conditions on the left endpoint
bc = DirichletBC(V,Constant('2.0'),q2b_boundary_left)

# Formulate the weak form 
u = TrialFunction(V)
v = TestFunction(V)
a = inner(q*grad(u),grad(v))*dx 
L = fa*v*dx + Constant('5.0')*q*v*ds

# Solve the problem
u = Function(V)
solve(a==L,u,bc)

# Make a plot
plot(u,interactive=True)