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

def quadrilateral ( ):

#*****************************************************************************80
#
## quadrilateral solves the Poisson problem using a quadrilateral mesh.
#
#  Discussion:
#
#    - div grad U = - 6        in Omega
#    U(dOmega) = Uexact(x,y)   on dOmega
#
#    Uexact = 1 + x^2 + 2 y^2
#    Omega = unit square [0,1]x[0,1]
#
#  Licensing:
#
#    This code is distributed under the MIT license.
#
#  Modified:
#
#    27 November 2018
#
#  Author:
#
#    John Burkardt
#
  import numpy as np
#
#  Create an 8x8 quadrilateral mesh on the unit square.
#  Documentation on the quadrilateral mesh option is scattered and
#  unreliable.
#
  mesh = UnitSquareMesh.create ( 64, 64, CellType.Type.quadrilateral )
# mesh = UnitSquareMesh ( 8, 8, CellType.Type_quadrilateral )
# mesh = UnitSquareMesh ( 8, 8, CellType.type.quadrilateral )
# mesh = UnitSquareMesh ( 8, 8, CellType.Type.quadrilateral )
# mesh = UnitQuadMesh ( 8, 8 )
#
#  Define the function space.
#
  V = FunctionSpace ( mesh, 'P', 1 )
#
#  Define the exact solution..
#
  uexact_expr = Expression ( '1 + x[0]*x[0] + 6*x[1]*x[1]', degree = 4 )
#
#  Define the boundary condition.
#
  def boundary ( x, on_boundary ):
    return on_boundary

  bc = DirichletBC ( V, uexact_expr, boundary )
#
#  Define the variational problem symbolically.
#
  u = TrialFunction ( V )
  v = TestFunction ( V )
  f = Constant ( -14.0 )
  a = dot ( grad ( u ), grad ( v ) ) * dx
  L = f * v * dx
#
#  Compute the solution.
#
  uh = Function ( V )
  solve ( a == L, uh, bc )
#
#  Plot the solution.
#  We will be able to plot the solution, and the mesh, using ParaView.
#
  filename = "quadrilateral.pvd"
  file = File ( filename )
  file << uh
  print ( '  Graphics saved as "%s"' % ( filename ) )
#
#  Compute maximum error at vertices.
#
  vertex_values_uexact = uexact_expr.compute_vertex_values ( mesh )
  vertex_values_uh = uh.compute_vertex_values ( mesh )
  error_max = np.max ( np.abs ( vertex_values_uexact - vertex_values_uh ) )
  print ( '  error_max =', error_max )
#
#  Terminate.
#
  return

def quadrilateral_test ( ):

#*****************************************************************************80
#
## quadrilateral_test tests quadrilateral.
#
#  Licensing:
#
#    This code is distributed under the MIT license.
#
#  Modified:
#
#    27 November 2018
#
#  Author:
#
#    John Burkardt
#
  import dolfin
  import platform
  import time

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

  print ( '' )
  print ( 'quadrilateral_test:' )
  print ( '  Python version: %s' % ( platform.python_version ( ) ) )
  print ( '  FENICS version %s'% ( dolfin.__version__ ) )
  print ( '  Poisson equation on the unit square with quad elements.' )

  quadrilateral ( )
#
#  Terminate.
#
  print ( '' )
  print ( 'quadrilateral_test:' )
  print ( '  Normal end of execution.' )
  print ( '' )
  print ( time.ctime ( time.time() ) )
  return

if ( __name__ == '__main__' ):

  quadrilateral_test ( )