# femth_bihar2dps_la.py  Galerkin FEM biharmonic 4th order  parallel
#                                                                   
# solve  uxxxx(x,y) + uyyyy(x,y) + 2*uxxyy(x,y) = f(x,y)
#
#  f(x,y) := 4.0*sin(x)*sin(y)
#
#  in the rectangle  xmin< x <xmax and  ymin< y <ymax  ,
#
#     subject to boundary conditions of the first
#     kind, Dirichlet
#
#  Dirichlet boundary conditions are computed using:
#  ub(x,y) := sin(x)*sin(y)
#
#  uses high order shape functions
#  uses high order Gauss-Legendre quadrature
#  used low degree of freedom
#  uses laphi.py  simeq.py  deriv.py  gauleg.py 
 

import math
import subprocess
import thread
import threading
from pybarrier import * # gets  barrier
from time import time
from laphi import *
from simeq import *
from deriv import *
from gauleg import *

global bar
global nx, ny, npx, npy, xg, yg, xx, wx, yy, wy

# PDE problem definition functions  f and ub 

def F(x, y):
  return 4.0*math.sin(x)*math.sin(y)

def ub(x, y):
  return math.sin(x)*math.sin(y)

def galk(x, y, i, j, ii, jj):
            # Galerkin k stiffness function for this problem       
            # solve  uxxxx(x,y) + 2*uxxyy(x,y) + uyyyy(x,y) = f(x,y) 
  return(( phipppp(x,j,nx-1,xg)*     phi(y,jj,ny-1,yg) +
        2.0* phipp(x,j,nx-1,xg)*   phipp(y,jj,ny-1,yg) +
               phi(x,j,nx-1,xg)* phipppp(y,jj,ny-1,yg) ) *
         (     phi(x,i,nx-1,xg)*     phi(y,ii,ny-1,yg) ) )

# end galk used by integration  

def galf(x, y, i, ii): # Galerkin f function for this problem 
  return F(x,y)*phi(x,i,nx-1,xg)*phi(y,ii,ny-1,yg)

# end galf used by integration  

def write_soln():
  print "writing femth_bihar2dps_la_py.dat  "
  fout = open("femth_bihar2dps_la_py.dat", "w")
  for i in range(nx):
    for ii in range(ny):
      fout.write(' '+str(xg[i])+'  '+str(yg[ii])+'  '+str(ug[i*ny+ii])+'\n')
 
    fout.write(" \n")

  fout.close
  print "finished writing femth_bihar2dps_la_py.dat"

# end write_soln 

def check_soln():
               # for all interior points 
  cxx = [[0 for j in range(nx)] for i in range(nx)]
  cxxxx = [[0 for j in range(nx)] for i in range(nx)]
  cyy = [[0 for j in range(ny)] for i in range(ny)]
  cyyyy = [[0 for j in range(ny)] for i in range(ny)]
  tcxxyy = [0 for j in range(ny)]
  tempx = [0 for j in range(nx)]
  tempy = [0 for j in range(ny)]

  cnt = 0
  for i in range(1,nx-1):
    tempx = rderiv(2,nx,i,hx)
    for k in range(0,nx):
      cxx[i][k] = tempx[k]
    tempx = rderiv(4,nx,i,hx)
    for k in range(0,nx):
      cxxxx[i][k] = tempx[k]

  for ii in range(1,ny-1):
    tempy = rderiv(2,ny,ii,hy)
    for k in range(0,ny):
      cyy[ii][k] = tempy[k]
    tempy = rderiv(4,ny,ii,hy)
    for k in range(0,ny):
      cyyyy[ii][k] = tempy[k]

  sumerr = 0.0
  sumsqerr = 0.0
  maxerr = 0.0

  for i in range(1,nx-1): # full equation 
    x = xmin + i*hx
    for j in range(1,ny-1):
      y = ymin + j*hy
      eval_deriv = 0.0 # no PDE u term 
        
      for k in range(0,nx):
	eval_deriv += cxxxx[i][k]*ug[k*ny+j]   # PDE  uxxxx term 
      for k in range(0,ny):
	eval_deriv += cyyyy[j][k]*ug[i*ny+k]    # PDE uyyyy term 

      #  + 2*uxxyy 
      for jj in range(0,ny): # walk jj in y, computing uxx 
	tcxxyy[jj] = 0.0 # numeric  uxx part of uxxyy 
        for k in range(0,nx):
	  tcxxyy[jj] += cxx[i][k]*ug[k*ny+jj]

      for k in range(0,ny):  # uyy part of uxxyy 
        eval_deriv += 2.0*cyy[j][k]*tcxxyy[k]  # PDE uxxyy term 

      err = eval_deriv-F(x,y)
      err = abs(err)
      sumerr += err
      sumsqerr += err*err
      maxerr = max(err,maxerr)
      cnt += 1  
    # end y 
  # end x 
  rmserr = math.sqrt(sumsqerr/cnt)
  avgerr = sumerr/cnt

  print "check_soln values"
  print "maxerr=",maxerr,", rmserr=",rmserr,", avgerr=",avgerr
    
# end check_soln  

class inner_loops(threading.Thread): # most work performed here
  def __init__ (self,myid):
    self.i = myid
    threading.Thread.__init__(self)
    print "inner_loops(",i,") init"

  def run(self):
    print "run(",self.i,") executing"
    # compute entries in stiffness matrix  
    # for i in range(1,nx-1): outer loop
    for j in range(0,nx):
      for ii in range(1,ny-1):
        for jj in range(0,ny):
          # compute integral for k(i,j)   
          val = 0.0
          for px in range(1,npx+1):
            for py in range(1,npy+1):
              val = val + wx[px]*wy[py]*galk(xx[px],yy[py],self.i,j,ii,jj)

          kg[(self.i*ny+ii)][(j*ny+jj)] = val

    # compute integral for f(i)   
    # for i in range(1,nx-1):  outer loop
    for ii in range(1,ny-1):
      val = 0.0
      for px in range(1,npx+1):
        for py in range(1,npy+1):
          val = val + wx[px]*wy[py]*galf(xx[px],yy[py],self.i,ii)

      fg[self.i*ny+ii] = val

    print "run(",self.i,") finished"
    bar.wait() # completion barrier in master

# end class inner_loops

nx = 7 #  9            # use  "global nx"  for classes
ny = 7 #  9
npx = 7 #  9
npy = 7 #  9

xmin = 0.0 # problem parameters 
xmax = 1.5 # 3.14159 
ymin = 0.0
ymax = 1.5 # 3.14159 

# nx, ny, xg, yg needed by galk or galf, available at call
xg = [0 for j in range(nx)] # X grid 
yg = [0 for j in range(ny)] # Y grid 
ug = [0 for j in range(nx*ny)] # solution 

nxy = nx*ny
# npx, npy  Gauss Legendre integration order 
xx = [0 for j in range(npx+1)]
wx = [0 for j in range(npx+1)]
yy = [0 for j in range(npy+1)]
wy = [0 for j in range(npy+1)]
kg = [[0 for j in range(nx*ny)] for i in range(nx*ny)]
fg = [0 for j in range(nx*ny)]
Ua = [0 for j in range(nx*ny)]

print "femth_bihar2dps_la.java running. "
now = time()
start = now
prev = now
print "solve  uxxxx(x,y) + 2*uxxyy(x,y) + uyyyy(x,y) = f(x,y) "
print "  "
print "f(x,y) := 4*sin(x)*sin(y)  "
print "  "
print "in the rectangle   xmin< x <xmax  and  ymin< y <ymax, "
print "  "
print "   subject to boundary conditions of the first kind, Dirichlet "
print "   ub(x,y) := sin(x)*sin(y)  "
print " using high order shape functions"
print " using high order quadrature"
print ((nx-2)*(ny-2))," Degrees of Freedom"
print "xmin=",xmin,", xmax=",xmax,", ymin=",ymin,", ymax=",ymax
print "nx=",nx,", ny=",ny
print "x grid and analytic solution at ymax"
hx = (xmax-xmin)/(nx-1)
for i in range(0,nx):
  xg[i] = xmin+i*hx
  print "i=",i,", xg[i]=",xg[i],", ub(xg[i],ymax)=",ub(xg[i],ymax)

print "  "
print "y grid and analytic solution at xmax"
hy = (ymax-ymin)/(ny-1)
for ii in range(0,ny):
  yg[ii] = ymin+ii*hy
  print "ii=",ii,", yg[ii]=",yg[ii],", ub(xmax,yg[ii])=",ub(xmax, yg[ii])

print "  "
# put solution in Ua vector  
for i in range(0,nx):
  x = xmin+i*hx
  for ii in range(0,ny):
    y = ymin+ii*hy
    Ua[i*ny+ii] = ub(x,y)

print "  "

# initialize k and f  
for i in range(0,nx):
  for j in range(0,nx):
    for ii in range(0,ny):
      for jj in range(0,ny):
        kg[(i*ny+ii)][(j*ny+jj)] = 0.0

for i in range(0,nx):
  for ii in range(0,ny):
    fg[i*ny+ii] = 0.0

# set boundary conditions, four sides  
for i in range(0,nx):
  x = xg[i]
  ii=0
  y = yg[ii]
  kg[(i*ny+ii)][(i*ny+ii)] = 1.0
  fg[i*ny+ii] = ub(x,y)
  ii=ny-1
  y = yg[ii]
  kg[(i*ny+ii)][(i*ny+ii)] = 1.0
  fg[i*ny+ii] = ub(x,y)

for ii in range(0,ny):
  y = yg[ii]
  i=0
  x = xg[i]
  kg[(i*ny+ii)][(i*ny+ii)] = 1.0
  fg[i*ny+ii] = ub(x,y)
  i=nx-1
  x = xg[i]
  kg[(i*ny+ii)][(i*ny+ii)] = 1.0
  fg[i*ny+ii] = ub(x,y)


print "calling gauleg xmin=",xmin,", xmax=",xmax,", npx=",npx
xx,wx = gaulegf(xmin, xmax, npx) # xx and wx set up for integrations  
print "xx[1]=",xx[1],", xx[2]=",xx[2],", wx[1]=",wx[1],", wx[2]=",wx[2]
print "calling gauleg ymin=",ymin,", ymax=",ymax,", npy=",npy
yy,wy = gaulegf(ymin, ymax, npy) # yy and wy set up for integrations 
print "yy[1]=",yy[1],", yy[2]=",yy[2],", wy[1]=",wy[1],", wy[2]=",wy[2]
i=1
ii=1
j=1
jj=1
val = galk(xx[2],yy[2],i,j,ii,jj)
print "galk(xx[2],yy[2],i,j,ii,jj)=",val
val = galf(xx[2],yy[2],i,ii)
print "galf(xx[2],yy[2],i,ii)=",val

now = time()
print "time to initialize = ",now-prev," seconds"
prev = now

# compute entries in stiffness matrix and RHS  
print "compute stiffness matrix  "
bar = barrier(nx-1) #  nx-2 threads plus main
for i in range(1,nx-1):
  inner_loops(i).start()

bar.wait() # all threads must finish to pass this barrier

print "k computed stiffness matrix"
print "f computed forcing function"
now = time()
print "time for stiffness and forcing = ",now-prev," seconds"
prev = now
ug = simeq(kg, fg)
now = time()
print "time for solve = ",now-prev," seconds"
prev = now

write_soln() # to file x y soln 

check_soln() # independent of nowing analytic solution 

avgerr = 0.0
maxerr = 0.0
print "ug computed Galerkin, Ua analytic, error  "
for i in range(0,nx):
  for ii in range(0,ny):
    err = abs(ug[i*ny+ii]-Ua[i*ny+ii])
    if err>maxerr:
      maxerr = err
    avgerr = avgerr + err

print "xmax=",xmax,", ymax=",ymax,", nx=",nx,", ny=",ny,
print ", npx=",npx,", npy=",npy
print "                maxerr=",maxerr,", avgerr=",avgerr/(nx*ny)
print " "

now = time()
print "total CPU time = ",now-start," seconds"
print "calling gnuplot"
stat = subprocess.call(["gnuplot","-persist","fem_bihar2dps_la_py.plot"])
print "returned stat= ", stat

# end femth_bihar2dps_la.py