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

def mitchell06 ( ):

#*****************************************************************************80
#
## mitchell06 sets up Mitchell test #6.
#
#  Boundary layer
#  -eps * Laplace(u) + 2 du/dx + du/dy = f on [-1,+1]x[-1,+1]
#  u = u0 on the boundary.
#  u0 = (1-exp(-(1-x)/eps)*(1-exp(-(1-y)/eps)*cos(pi(x+y))
#  f = -eps*(u0'')+2 du0/dx+du0/dy
#
#  Discussion:
#
#    Suggested parameter values:
#    * eps = 0.1
#    * eps = 0.001
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license.
#
#  Modified:
#
#    30 August 2014
#
#  Author:
#
#    John Burkardt
#
#  Reference:
#
#    William Mitchell,
#    A collection of 2D elliptic problems for testing adaptive
#    grid refinement algorithms,
#    Applied Mathematics and Computation,
#    Volume 220, 1 September 2013, pages 350-364.
#
  import matplotlib.pyplot as plt
#
#  Set parameters.
#
  eps = 0.1
#
#  Define a mesh of triangles on [-1,+1]x[-1,+1]
#
  mesh = RectangleMesh ( -1.0, -1.0, +1.0, +1.0, 4, 4 )
#
#  Plot the mesh.
#
  plot ( mesh, title = 'Mitchell Test 06 Mesh' )
  filename = 'mitchell06_mesh.png'
  plt.savefig ( filename )
  print ( '  Graphics saved as "%s"' % ( filename ) )
  plt.close ( )
#
#  Define the function space.
#
  V = FunctionSpace ( mesh, "Lagrange", 1 )
#
#  Define the exact solution.
#
  u0 = Expression ( "( 1.0 - exp ( - ( 1.0 - x[0] ) / eps ) ) * " 
                    "( 1.0 - exp ( - ( 1.0 - x[1] ) / eps ) ) * " 
                    "cos ( pi * ( x[0] + x[1] ) )", degree = 10 )
#
#  Apply the boundary condition to points which are on the boundary of the mesh.
#
  def u0_boundary ( x, on_boundary ):
    return on_boundary
#
#  The value to be used at Dirichlet boundary points is U0.
#
  bc = DirichletBC ( V, u0, u0_boundary )
#
#  Define test and trial functions.
#
  u = TrialFunction ( V )
  v = TestFunction ( V )
#
#  Define the right hand side.
#
  f = Expression \
    ( "exp ( - 2.0 / eps ) / pow ( eps, 2 ) * (" 
      " - exp ( ( 1.0 - x[1] ) / eps )" 
      " * ( -2.0 * exp ( 1.0 / eps ) * pow ( pi * eps, 2 )"
      " + exp ( x[0] / eps ) * ( 1.0 + 2.0 * pow ( pi * eps, 2 ) ) ) "
      " * cos ( pi * ( x[0] + x[1] ) ) "
      " - ( 3.0 * exp ( 2.0 / eps ) - exp ( ( 1.0 - x[0] ) / eps )"
      "  - exp ( ( 1.0 - x[1] ) / eps ) - exp ( ( x[0] + x[1] ) / eps ) )"
      "  * eps * pi * sin ( pi * ( x[0] + x[1] ) ) )",
      eps = eps,
      pi = pi, degree = 10
    )
#
#  Define the variational problem.
#  How can I BREAK UP THIS LONG LINE?
#  How do I express 'DUDX'?
#
  a = ( eps * inner ( nabla_grad ( u ), nabla_grad ( v ) ) + 2.0 * inner ( DUDX, v ) + inner ( DUDY, v ) ) * dx
  L = f * v * dx
#
#  Compute the solution.
#
  u = Function ( V )
  solve ( a == L, u, bc )
#
#  Plot the solution.
#
  plot ( u )

  return

def mitchell06_test ( ):

#*****************************************************************************80
#
## mitchell06_test tests mitchell06.
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license.
#
#  Modified:
#
#    22 October 2018
#
#  Author:
#
#    John Burkardt
#

#
#  Report level = only warnings or higher.
#
  level = 30
  set_log_level ( level )

  print ( '' )
  print ( 'mitchell06_test:' )
  print ( '  FENICS/Python version' )
  print ( '  Mitchell test problem #6.' )

  mitchell06 ( )
#
#  Terminate.
#
  print ( '' )
  print ( 'mitchell06_test:' )
  print ( '  Normal end of execution.' )
  return

if ( __name__ == '__main__' ):

  mitchell06_test ( )