import sys, time
import numpy as np
seed = 1435
np.random.seed(seed)

# import problem classes
import Exeter_CFD_Problems as TestProblems
verbose = False

sys.argv = sys.argv[:1]
print('Demonstration of the PitzDialy test problem.')
# set up directories.
settings = {
    'source_case': 'Exeter_CFD_Problems/data/PitzDaily/case_fine/',
    'case_path': 'Exeter_CFD_Problems/data/PitzDaily/case_single/',
    'boundary_files': ['Exeter_CFD_Problems/data/PitzDaily/boundary.csv'],
    'fixed_points_files': ['Exeter_CFD_Problems/data/PitzDaily/fixed.csv']
}
# instantiate the problem object
prob = TestProblems.PitzDaily(settings)
# get the lower and upper bounds
lb, ub = prob.get_decision_boundary()
# generate random solutions satisfying the lower and upper bounds.
x = np.array([ 1.54630729e-01, -4.99600344e-02, 1.38160660e-01, -4.99666603e-02,
2.28714652e-04, -3.24888531e-04, 2.74147357e-01, -1.88668377e-02,
6.62354130e-04, 1.86793810e-04])
#rand_x = []
#for i in range(x.shape[0]):
#    if prob.constraint(x[i]): # check to see if the random solution is valid
#    rand_x.append(x[i])
# evaluate for a solution
print('Number of control points: ', prob.n_control)
print('Decision vector: ', x)
print('Running simulation ...')
start = time.time()
res = prob.evaluate(x, verbose=verbose)
print('Objective function value:', res)
print('Time taken:', time.time()-start, ' seconds.')

print(max(lb),min(lb))

print(max(ub),min(ub))

x

! ls

print(case.fields)

case.xlim[0]

case.ylim[0]

case.xlim[1]

max(case.cellCentres[0])

case.cellCentres

case['U']

case["magU"]

plot=plot_field(case, "magU", cmap='rainbow',colorbar=False)
plt.xlim([-0.025,0.2])
plt.ylim([-0.03,0.03])
plt.xlabel(r'$x \ (m)$')
plt.ylabel(r'$y \ (m)$')
bar=add_colorbar(plot, aspect=10, padFraction=1)
bar.set_label(r'$Velocity \ Magnitude \ (m/s)$')

plot2=plot_vectors(case,'U', colorField='magU',cmap='rainbow', normalize=True, \
             sampleByPlane=True,planeResolution=[20, 80], width=0.0015)
plt.xlim([-0.025,0.2])
plt.ylim([-0.03,0.03])
plt.xlabel(r'$x \ (m)$')
plt.ylabel(r'$y \ (m)$')
bar2=add_colorbar(plot2, aspect=10, padFraction=1)
bar2.set_label(r'$Velocity \ Magnitude \ (m/s)$')

plot3=plot_streamlines(case, 'U',colorField='magU',cmap='rainbow',planeResolution=[100,100],density=2)
plt.xlim([-0.025,0.2])
plt.ylim([-0.03,0.03])
plt.xlabel(r'$x \ (m)$')
plt.ylabel(r'$y \ (m)$')
bar3=add_colorbar(plot3.lines, aspect=10, padFraction=1)
bar3.set_label(r'$Velocity \ Magnitude \ (m/s)$')

plot4=plot_streamlines(case2, 'U',colorField='magU',cmap='rainbow',planeResolution=[90,90],density=2)
plt.xlim([-0.025,0.2])
plt.ylim([-0.03,0.03])
plt.xlabel(r'$x \ (m)$')
plt.ylabel(r'$y \ (m)$')
bar4=add_colorbar(plot4.lines, aspect=10, padFraction=1)
bar4.set_label(r'$Velocity \ Magnitude \ (m/s)$')

# Partial re-run of the case to access values for plotting
import sys, time
import numpy as np
seed = 1435
np.random.seed(seed)

# import problem classes
import Exeter_CFD_Problems as TestProblems
verbose = False

sys.argv = sys.argv[:1]
print('Demonstration of the PitzDialy test problem.')
# set up directories.
settings = {
    'source_case': 'Exeter_CFD_Problems/data/PitzDaily/case_fine/',
    'case_path': 'Exeter_CFD_Problems/data/PitzDaily/case_single/',
    'boundary_files': ['Exeter_CFD_Problems/data/PitzDaily/boundary.csv'],
    'fixed_points_files': ['Exeter_CFD_Problems/data/PitzDaily/fixed.csv']
}
# instantiate the problem object
prob = TestProblems.PitzDaily(settings)
# get the lower and upper bounds
lb, ub = prob.get_decision_boundary()
# generate random solutions satisfying the lower and upper bounds.
x = np.random.random((1000, lb.shape[0])) * (ub - lb) + lb
rand_x = []
for i in range(x.shape[0]):
    if prob.constraint(x[i]): # check to see if the random solution is valid
        rand_x.append(x[i])

boundaries=plot_boundaries(case2)
plt.xlim([-0.025,0.3])
plt.ylim([-0.055,0.03])
plt.xlabel(r'$x \ (m)$')
plt.ylabel(r'$y \ (m)$')

# decision vector
for i in range(0,10,2):
    plt.plot(rand_x[0][i],rand_x[0][i+1], 'r^')

# boundary points
pts = [np.loadtxt(filename, delimiter=',') for filename in prob.domain_files]

x_boundary = []
y_boundary = []

for i in range(0,11,1):
    x_boundary.append(pts[0][i][0])
    y_boundary.append(pts[0][i][1])

plt.plot(x_boundary,y_boundary,'bo--')

#fixed_points
fixed_points = [list(np.loadtxt(filename, delimiter=',').astype(int))\
                            for filename in prob.fixed_points_files]

fixed=[fixed_points[0][0][0],fixed_points[0][0][1],fixed_points[0][1][0],fixed_points[0][1][1]]

x_fixed = []
y_fixed = []

for j in fixed:
    x_fixed.append(pts[0][j][0])
    y_fixed.append(pts[0][j][1])

plt.plot(x_fixed,y_fixed,'ks')

x_curve = []
y_curve = []

#catmull_clark curve
shape=prob.convert_decision_to_shape(rand_x[0]) # convert random array to ordered list
for control_points in shape:
    data_points = TestProblems.data.support.catmull(control_points)
    for i in range(prob.niter):
        data_points = TestProblems.data.support.catmull(data_points)

for k in range(0,512):
    x_curve.append(data_points[k][0])
    y_curve.append(data_points[k][1])

plt.plot(x_curve,y_curve,'g-')

len(pts[0])

pts #these are the boundary points [x, y] there are 11 of them

fixed_points # these are the fixed points, the first two points, and points 4 and 5. They are grouped like
             # this to control the gradient at these points

! foamToVTK -case Exeter_CFD_Problems/data/PitzDaily/case_non_optimised > log3 2>&1
grid = \
'Exeter_CFD_Problems/data/PitzDaily/case_non_optimised/VTK/case_non_optimised_1634.vtk'
plotPitzDailyGrid.plot(grid)

! foamToVTK -case Exeter_CFD_Problems/data/PitzDaily/case_single > log4 2>&1
grid_2 = \
'Exeter_CFD_Problems/data/PitzDaily/case_single/VTK/case_single_0.vtk'
plotPitzDailyGrid.plot(grid_2)