from firedrake import *
import numpy as np

mesh = UnitSquareMesh(100,100)
Q = FunctionSpace(mesh,'DG',0)

file_prefix = 'perm600'
class createData(Expression):
  def __init__(self):
    self.s = np.loadtxt(file_prefix+'_s')
    self.xs = np.loadtxt(file_prefix+'_xs')
    self.ys = np.loadtxt(file_prefix+'_ys')
  def eval(self, val, x):
    numerator = 0
    denominator = 0
    for i in range(0,600):
      numerator += self.s[i]/((self.xs[i]-x[0])**2+(self.ys[i]-x[1])**2)
      denominator += 1/((self.xs[i]-x[0])**2+(self.ys[i]-x[1])**2)
    ssum = numerator/denominator
    val[0] = 1e-12*exp(ssum)
  def value_shape(self):
    return ()

u = Function(Q)
u.interpolate(createData())