#! /usr/bin/env python
#
def fem1d_bvp_linear_test07 ( ):

#*****************************************************************************80
#
## FEM1D_BVP_LINEAR_TEST07 does an error analysis.
#
#  Location:
#
#    http://people.sc.fsu.edu/~jburkardt/py_src/fem1d_bvp_linear/fem1d_bvp_linear_test07.py
#
#  Discussion:
#
#    Use A7, C7, F7, EXACT7, EXACT_UX7.
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license.
#
#  Modified:
#
#    09 January 2015
#
#  Author:
#
#    John Burkardt
#
#  Reference:
#
#    Eric Becker, Graham Carey, John Oden,
#    Finite Elements, An Introduction, Volume I,
#    Prentice-Hall, 1981, page 123-124,
#    ISBN: 0133170578,
#    LC: TA347.F5.B4.
#
  import numpy as np
  import platform
  from fem1d_bvp_linear import fem1d_bvp_linear
  from h1s_error_linear import h1s_error_linear
  from l1_error import l1_error
  from l2_error_linear import l2_error_linear
  from max_error_linear import max_error_linear

  print ( '' )
  print ( 'FEM1D_BVP_LINEAR_TEST07' )
  print ( '  Python version: %s' % ( platform.python_version ( ) ) )
  print ( '  Becker/Carey/Oden example.' )
  print ( '  Solve -( A(x) U\'(x) )\' + C(x) U(x) = F(x)' )
  print ( '  for 0 < x < 1, with U(0) = U(1) = 0.' )
  print ( '' )
  print ( '  Compute L2 norm and seminorm of error for various N.' )
  print ( '' )
  print ( '     N        L1 error        L2 error      Seminorm error  Maxnorm error' )
  print ( '' )

  n = 11

  for i in range ( 0, 5 ):
#
#  Geometry definitions.
#
    x_first = 0.0
    x_last = 1.0
    x = np.linspace ( x_first, x_last, n )

    u = fem1d_bvp_linear ( n, a7, c7, f7, x )

    e1 = l1_error ( n, x, u, exact7 )
    e2 = l2_error_linear ( n, x, u, exact7 )
    h1s = h1s_error_linear ( n, x, u, exactp7 )
    mx = max_error_linear ( n, x, u, exact7 )

    print ( '  %4d  %14g  %14g  %14g  %14g' % ( n, e1, e2, h1s, mx ) )

    n = 2 * ( n - 1 ) + 1
#
#  Terminate.
#
  print ( '' )
  print ( 'FEM1D_BVP_LINEAR_TEST07' )
  print ( '  Normal end of execution.' )
  return

def a7 ( x ):
  alpha = 30.0
  x0 = 1.0 / 3.0
  value = 1.0 / alpha + alpha * ( x - x0 ) ** 2
  return value

def c7 ( x ):
  value = 0.0
  return value

def exact7 ( x ):
  from math import atan
  alpha = 30.0
  x0 = 1.0 / 3.0
  value = ( 1.0 - x ) * ( atan ( alpha * ( x - x0 ) ) + atan ( alpha * x0 ) )
  return value

def exactp7 ( x ):
  from math import atan
  alpha = 30.0
  x0 = 1.0 / 3.0
  value = - atan ( alpha * ( x - x0 ) ) - atan ( alpha * x0 ) \
    + ( 1.0 - x ) * alpha / ( 1.0 + alpha * alpha * ( x - x0 ) ** 2 )
  return value

def f7 ( x ):
  from math import atan
  alpha = 30.0
  x0 = 1.0 / 3.0
  value = 2.0 * ( 1.0 + alpha * ( x - x0 ) * \
    ( atan ( alpha * ( x - x0 ) ) + atan ( alpha * x0 ) ) )
  return value

if ( __name__ == '__main__' ):
  from timestamp import timestamp
  timestamp ( )
  fem1d_bvp_linear_test07 ( )
  timestamp ( )