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

def poisson_mixed ( my_degree ):

#*****************************************************************************80
#
## poisson_mixed solves a Poisson equation using a mixed method.
#
#  Discussion:
#
#    This code is an adaptation of the standard fenics file demo_mixed-poisson.py.
#
#        sigma - grad u = 0    in Omega
#    div sigma          = - f  in Omega
#
#    u           = 0           on Gamma_D
#    sigma dot n = sin(5*x)    on Gamma_N
#
#  Licensing:
#
#    This code is distributed under the MIT license.
#
#  Modified:
#
#    11 January 2019
#
#  Author:
#
#    John Burkardt
#
#  Parameters:
#
#    Input, integer MY_DEGREE, specifies the approximation degree.
#    BDM, for sigma, uses MY_DEGREE+1,
#    DG, for u, uses MY_DEGREE.
#
  import matplotlib.pyplot as plt
#
#  Create the mesh.
#
  mesh = UnitSquareMesh ( 32, 32 )
#
#  Define finite elements spaces and build mixed space
#
  print ( '  BDM degree = %d, DG degree = %d' % ( my_degree + 1, my_degree ) )

  BDM = FiniteElement ( "BDM", mesh.ufl_cell(), my_degree + 1 )
  DG  = FiniteElement ( "DG", mesh.ufl_cell(), my_degree )
#
#  Build the mixed space.
#
  W = FunctionSpace ( mesh, BDM * DG )
#
#  Specify the trial functions (the unknowns) and the test functions.
#
  ( sigma, u ) = TrialFunctions ( W )
  ( tau, v ) = TestFunctions ( W )
#
#  Define the source function
#
  f = Expression ( "10*exp(-(pow(x[0] - 0.5, 2) + pow(x[1] - 0.5, 2)) / 0.02)", degree = 2 )
#
#  Define the variational form
#
  a = ( dot ( sigma, tau ) + div ( tau ) * u + div ( sigma ) * v ) * dx
  L = - f * v * dx
#
#  Define function G such that G \cdot n = g
#
  class BoundarySource ( UserExpression ):
    def __init__(self, mesh, **kwargs):
      self.mesh = mesh
      super().__init__(**kwargs)
    def eval_cell(self, values, x, ufc_cell):
      cell = Cell(self.mesh, ufc_cell.index)
      n = cell.normal(ufc_cell.local_facet)
      g = sin ( 5 * x[0] )
      values[0] = g*n[0]
      values[1] = g*n[1]
    def value_shape(self):
      return (2,)

  G = BoundarySource ( mesh, degree = 2 )
#
#  Define essential boundary
#
  def boundary ( x ):
    return x[1] < DOLFIN_EPS or x[1] > 1.0 - DOLFIN_EPS

  bc = DirichletBC ( W.sub(0), G, boundary )
#
#  Solve.
#
  w = Function ( W )
  solve ( a == L, w, bc )
#
#  Split the solution into components sigma and u.
#
  ( sigma, u ) = w.split()
#
#  Plot sigma and u.
#
  plot ( sigma, mesh = mesh, title = 'sigma' )
  filename = 'sigma.png'
  plt.savefig ( filename )
  print ( '  Graphics saved as "%s"' % ( filename ) )
  plt.close ( )

  plot ( u, mesh = mesh, title = 'u' )
  filename = 'u.png'
  plt.savefig ( filename )
  print ( '  Graphics saved as "%s"' % ( filename ) )
  plt.close ( )
#
#  Terminate.
#
  return

def poisson_mixed_test ( ):

#*****************************************************************************80
#
## poisson_mixed_test tests poisson_mixed.
#
#  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 ( 'poisson_mixed_test:' )
  print ( '  Python version: %s' % ( platform.python_version ( ) ) )
  print ( '  FENICS version %s'% ( dolfin.__version__ ) )
  print ( '  Mixed method for Poisson problem.' )
#
#  Report level = only warnings or higher.
#
  level = 30
  set_log_level ( level )

  my_degree = 1
  poisson_mixed ( my_degree )
#
#  Terminate.
#
  print ( '' )
  print ( 'poisson_mixed_test:' );
  print ( '  Normal end of execution.' )
  print ( '' )
  print ( time.ctime ( time.time() ) )

if ( __name__ == '__main__' ):

  poisson_mixed_test ( )