from firedrake import *
import numpy as np

mesh = UnitSquareMesh(100,100,quadrilateral=True)
Q = FunctionSpace(mesh,'DG',0)
xQ = VectorFunctionSpace(mesh,'DG',0)

file_prefix = 'perm600'
datas = np.loadtxt(file_prefix+'_s')
dataxs = np.loadtxt(file_prefix+'_xs')
datays = np.loadtxt(file_prefix+'_ys')
u = Function(Q)

x = SpatialCoordinate(Q.mesh())
x_DG = project(x,xQ)
usize = u.dat.data.size
for j in range(0,usize):
  numerator = 0
  denominator = 0
  for i in range(0,600):
    numerator += datas[i]/(4*np.power(dataxs[i]-x_DG.dat.data_ro[j,0],2) +np.power(datays[i]-x_DG.dat.data_ro[j,1],2))
    denominator += 1/(4*np.power(dataxs[i]-x_DG.dat.data_ro[j,0],2) +np.power(datays[i]-x_DG.dat.data_ro[j,1],2))
  value = 1e-12*np.exp(numerator/denominator)
  u.dat.data[j] = value
    
outfile = File('test.pvd')
outfile.write(u)