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

def dg_advection_diffusion ( my_grid, my_degree ):

#*****************************************************************************80
#
## dg_advection_diffusion : discontinuous Galerkin method, advection-diffusion.
#
#  Discussion:
#
#    u grad phi - kappa del phi = f
#
#    Dirichlet boundary on right side with value sin(5*pi*y).
#    Convecting velocity (-1, -0.4).
#
#  Licensing:
#
#    This code is distributed under the MIT license.
#
#  Modified:
#
#    11 January 2019
#
#  Author:
#
#    John Burkardt
#
#  Parameters:
#
  import matplotlib.pyplot as plt
#
#  Report input.
#
  print ( '' )
  print ( '  Case of %d x %d mesh' % ( my_grid, my_grid ) )
  print ( '  Approximation degree %d:' % ( my_degree ) )
#
#  Mesh the unit square.
#
  my_mesh = UnitSquareMesh ( my_grid, my_grid, 'right/left' )
#
#  Set function spaces.
#
  V_dg = FunctionSpace ( my_mesh, 'DG',  my_degree )
  V_cg = FunctionSpace ( my_mesh, 'CG',  my_degree )
  V_u  = VectorFunctionSpace ( my_mesh, 'CG', my_degree + 1 )
#
#  Create velocity function.
#
  u = interpolate ( Constant ( ( "-1.0", "-0.4" ) ), V_u )
#
#  Trial and test functions.
#
  phi = TrialFunction ( V_dg )
  v   = TestFunction ( V_dg )
#
#  Diffusivity.
#
  kappa = Constant ( 0.0 )
#
#  Source term.
#
  f = Constant ( 0.0 )
#
#  Penalty term.
#
  alpha = Constant ( 5.0 )
#
#  Mesh-related functions.
#
  n = FacetNormal ( my_mesh )
  h = CellDiameter ( my_mesh )
  h_avg = ( h('+') + h('-') ) / 2.0
#
#  (dot(v,n) + |dot(v,n)| ) / 2.0
#
  un = ( dot ( u, n ) + abs ( dot ( u, n ) ) ) / 2.0
#
#  Bilinear form.
#
  a_int = dot ( grad ( v ), kappa * grad ( phi ) - u * phi ) * dx

  a_fac = kappa('+') * ( alpha('+') / h('+') ) * dot ( jump ( v, n ), jump ( phi, n ) ) * dS \
        - kappa('+') * dot ( avg ( grad ( v ) ), jump ( phi, n ) ) * dS \
        - kappa('+') * dot ( jump ( v, n ), avg ( grad ( phi ) ) ) * dS

  a_vel = dot ( jump ( v ), un('+') * phi('+') - un('-') * phi('-') ) * dS \
        + dot ( v, un * phi ) * ds

  a = a_int + a_fac + a_vel
#
#  Linear form.
#
  L = v * f * dx
#
#  Set up boundary condition, applying strong BCs.
#
  class DirichletBoundary ( SubDomain ):
    def inside ( self, x, on_boundary ):
      return abs ( x[0] - 1.0 ) < DOLFIN_EPS and on_boundary

  g = Expression ( "sin ( pi * 5.0 * x[1] )", degree = 2 )
  bc = DirichletBC ( V_dg, g, DirichletBoundary(), "geometric" )
#
#  Solution function.
#
  phi_h = Function ( V_dg )
#
#  Assemble and apply boundary conditions.
#
  A = assemble ( a )
  b = assemble ( L )
  bc.apply ( A, b )
#
#  Solve system.
#
  solve ( A, phi_h.vector(), b )
#
#  Project solution to a continuous function space.
#
  up = project ( phi_h, V = V_cg )
#
#  Plot the solution.
#
  plot ( up, title = 'up' )
  filename = 'up.png'
  plt.savefig ( filename )
  print ( '  Graphics saved as "%s"' % ( filename ) )
  plt.close ( )
#
#  Terminate.
#
  return

def dg_advection_diffusion_test ( ):

#*****************************************************************************80
#
## dg_advection_diffusion_test tests dg_advection_diffusion.
#
#  Licensing:
#
#    This code is distributed under the MIT license.
#
#  Modified:
#
#    11 January 2019
#
#  Author:
#
#    John Burkardt
#
  import dolfin
  import platform
  import time

  print ( time.ctime ( time.time() ) )
  print ( '' )
  print ( 'dg_advection_diffusion_test:' )
  print ( '  Python version: %s' % ( platform.python_version ( ) ) )
  print ( '  FENICS version %s'% ( dolfin.__version__ ) )
  print ( '  Discontinuous Galerkin method for advection/diffusion equation.' )
#
#  Report level = only warnings or higher.
#
  level = 30
  set_log_level ( level )

  my_grid = 64
  my_degree = 1
  dg_advection_diffusion ( my_grid, my_degree )
#
#  Terminate.
#
  print ( '' )
  print ( 'dg_advection_diffusiontest:' );
  print ( '  Normal end of execution.' )
  print ( '' )
  print ( time.ctime ( time.time() ) )

if ( __name__ == '__main__' ):

  dg_advection_diffusion_test ( )