#! /usr/bin/env python
#
from dolfin import *

def bvp_01 ( e_num ):

#*****************************************************************************80
#
## BVP_01, 1D test problem for FENICS.
#
#  Discussion:
#
#    -u'' + u = x, 0 < x < 1
#    u(0) = 0, u(1) = 0
#
#    Exact solution is u(x) = x - sinh(x)/sinh(1)
#
#  Location:
#
#    http://people.sc.fsu.edu/~jburkardt/examples/pbs/bvp_01.py
#    http://people.sc.fsu.edu/~jburkardt/examples/pbs/bvp_01_mesh.xml
#    http://people.sc.fsu.edu/~jburkardt/examples/pbs/bvp_01_u.xml
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license.
#
#  Modified:
#
#    19 February 2016
#
#  Author:
#
#    John Burkardt
#
#  Parameters:
#
#    Input, integer E_NUM, the number of elements.
#
  print ""
  print "BVP_01:"
  print "  Solve  -u'' + u = x, 0 < x < 1"
  print "  u(0) = 0, u(1) = 0"
  print "    Exact solution is u(x) = x - sinh(x)/sinh(1)"
#
#  Create a mesh on the unit interval.
#
  mesh = UnitIntervalMesh ( e_num )
#
#  Define the function space to be of Lagrange type
#  using piecewise linear basis functions.
#
  V = FunctionSpace ( mesh, "Lagrange", 1 )
#
#  For convenience, we define the exact solution.
#
  exact = Expression ( "x[0] - sinh ( x[0] ) / sinh ( 1.0 )" )
#
#  Define U0_BOUNDARY ( ) to indicate boundary points.
#
  def u0_boundary ( x, on_boundary ):
    return on_boundary
#
#  The boundary conditions, stored as "bc", are Dirichlet boundary conditions,
#  defined using the function space V, with value U0, for points for which 
#  U0_BOUNDARY is true.
#
  bc = DirichletBC ( V, exact, u0_boundary )
#
#  The "TEST" function is the one that multiplies the PDE.
#  The "TRIAL" function is the one that is used to build the solution.
#
  psii = TestFunction ( V )
  psij = TrialFunction ( V )
#
#  Define the system matrix, A = integral ( psii' * psij' + psii * psij ) dx
#
  A = ( inner ( nabla_grad ( psii ), nabla_grad ( psij ) ) + psii * psij ) * dx
#
#  Define the function on the right hand side.
#
  f = Expression ( "x[0]" )
#
#  Define the linear system right hand side, RHS = integral ( psii * f ) dx
#
  RHS = psii * f * dx
#
#  Specify that the solution should be a linear combination of elements of V.
#
  u = Function ( V )
#
#  Solve the problem for U, with the given boundary conditions.
#
  solve ( A == RHS, u, bc )
#
#  Dump the mesh to an XML file.
#
  mesh_file = File ( "bvp_01_mesh.xml" )
  mesh_file << mesh
#
#  Dump the solution to an XML file.
#
  solution_file = File ( "bvp_01_u.xml" )
  solution_file << u
#
#  Plot the solution.
#  We can't plot to the screen if running a batch job.
#
# plot ( u,  )
#
#  The interactive command pauses the program, which for us means
#  that the plot stays on the screen long enough to admire.
#
# interactive ( )
#
#  Terminate.
#
  print ""
  print "BVP_01:"
  print "  Normal end of execution."
  print ""
  return

if ( __name__ == '__main__' ):
  from timestamp import timestamp
  timestamp ( )
  bvp_01 ( 8 )
  timestamp ( )