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

def interpolant ( ):

#*****************************************************************************80
#
## interpolant shows how to read in a mesh and an interpolant function.
#
#  Discussion:
#
#    The solution U of a finite element linear system is a kind of function,
#    which we can evaluate at a point X using an expression like "u(x)".
#
#    It is often convenient to define other functions for use in a finite 
#    element calculation.  Such a function is typically actually a "finite 
#    element interpolant", defined by giving its value at every node on the 
#    finite element mesh.
#
#    When a piecewise linear mesh is being used, the definition of the mesh 
#    (the vertices) specifies the location of all the nodes (or degrees
#    of freedom) so we can create an interpolant if we can list a function value
#    to associate with each vertex listed in the mesh definition.
#
#    In this example, a 2D triangular mesh has been defined in 
#    inteporlant_mesh.xml.
#
#    For each node in that mesh, a function value is supplied in the text file
#    interpolant_values.txt.
#
#    This script shows how to read the mesh, and then create a function q 
#    by reading in the function values.  The resulting function is plotted, to
#    show that it works just like a finite element solution function.
#
#    While the data can be completely arbitrary, in this case, the node
#    data was generated by sampling the formula
#      f(x,y) = cos(pi*x) * sin(2*pi*y)
#
#    Therefore, the interpolant function q(x,y) can be compared to f(x,y)
#    over the region.
#
#  Licensing:
#
#    This code is distributed under the MIT license.
#
#  Modified:
#
#    23 October 2018
#
#  Author:
#
#    Hans-Werner van Wyk
#
  from math import pi
  import numpy as np
  import matplotlib.pyplot as plt
#
#  Read an XML file describing a mesh based on a 30x30 grid of points
#  over the unit square.
# 
  mesh = Mesh ( 'interpolant_mesh.xml' ) 
#
#  Plot the mesh.
#
  plot ( mesh, title = 'User-defined mesh on unit square' )
  filename = 'interpolant_mesh.png'
  plt.savefig ( filename )
  print ( '  Graphics saved as "%s"' % ( filename ) )
  plt.close ( )
#
#  Define a function space based on the mesh.
#
  V = FunctionSpace ( mesh, "CG", 1 )
# 
#  Open a text file, containing a scalar value at each of the 900 nodes.
#
  file_vals = open ( 'interpolant_values.txt', 'r' )
#
#  Read each line of the file into a list.
#
  list_vals = [ float ( line ) for line in file_vals.readlines() ]
#
#  Convert the list to a numeric array.
#
  array_vals = np.array ( list_vals )
#
#  Let Q be a function in the function space V.
#
  q = Function ( V )
#
#  Get the values of Q from the numeric array.
#  FENICS uses a different ordering for the data than the mesh does.
#
  d2v = dof_to_vertex_map ( V )
  q.vector()[:] = array_vals[d2v]
#
#  Plot the Q function.
#
  plot ( q, title = 'Q(X,Y) = interpolant to F(X,Y) data.' )
  filename = 'interpolant_values.png'
  plt.savefig ( filename )
  print ( '  Graphics saved as "%s"' % ( filename ) )
  plt.close ( )
#
#  Define the function F(X,Y) which generated this data.
#
  f = Expression ( "cos ( pi * x[0] ) * sin ( 2 * pi * x[1] )", \
    pi = pi, degree = 10 );
#
#  Evaluate F on the mesh.
#
  fmesh = interpolate ( f, V )
#
#  Plot F.
#
  plot ( fmesh, title = 'F(X,Y) = cos(pi x) sin ( 2 pi y )' )
  filename = 'formula_values.png'
  plt.savefig ( filename )
  print ( '  Graphics saved as "%s"' % ( filename ) )
  plt.close ( )
#
#  Terminate.
#
  return

def interpolant_test ( ):

#*****************************************************************************80
#
## interpolant_test tests interpolant.
#
#  Licensing:
#
#    This code is distributed under the MIT license.
#
#  Modified:
#
#    23 October 2018
#
#  Author:
#
#    John Burkardt
#
  import dolfin
  import platform
#
#  Report level = only warnings or higher.
#
  level = 30
  set_log_level ( level )

  print ( '' )
  print ( 'interpolant_test:' )
  print ( '  Python version: %s' % ( platform.python_version ( ) ) )
  print ( '  FENICS version %s'% ( dolfin.__version__ ) )
  print ( '  Define a function by reading mesh and value files.' )

  interpolant ( )
#
#  Terminate.
#
  print ( '' )
  print ( 'interpolant_test:' )
  print ( '  Normal end of execution.' )
  return

if ( __name__ == '__main__' ):
  import time
  print ( time.ctime ( time.time() ) )
  interpolant_test ( )
  print ( time.ctime ( time.time() ) )