using JFVM

Nx = 10 # number of cells in x direction
Lx = 5.0 # [m] domain length in x direction
m_uni = createMesh1D(Nx, Lx)

m_uni.facecenters.x

x = [0.0, 0.1, 0.3, 0.35, 1.0, 3.0, 3.4, 4.5, 5.0]
m_nonuni = createMesh1D(x)

Nx = 3 # number of cells in x-direction
Ny = 4 # number of cells in y-direction
Nz = 5 # number of cells in z-direction
Nr = 6 # number of cells in r-direction [radial]
Ntheta = 7 # number of cells in \theta-direction [2d radial, 3d cylindrical and spherical]

Lx = 1.0 # domain length in x-direction
Ly = 2.0 # domain length in y-direction
Lz = 3.0 # domain length in z-direction
Lr = 4.0 # domain length in r-direction
Ltheta = 2*pi # domain length in \theta-direction [in radian]

x=[0.0 , 0.4, 0.6, 1.0] # cell face positions in the x-direction
y=[0.0 , 0.4, 0.6, 1.0, 2.0] # cell face positions in the y-direction
z=[0.0 , 0.4, 0.6, 1.0, 2.0, 3.0] # cell face positions in the z-direction
r=[0.0 , 0.4, 0.6, 1.0, 2.0, 3.0, 4.0] # cell face positions in the r-direction
theta=5.0./[0.0 , 0.4, 0.6, 1.0, 2.0, 3.0, 4.0, 5.0]*Ltheta # cell face positions in the \theta-direction

# 1D cartesian uniform mesh:
m_1d_cartesian_uniform = createMesh1D(Nx, Lx)
# 1D cartesian nonuniform mesh
m_1d_cartesian_nonuniform = createMesh1D(x)
# 1D radial uniform mesh
m_1d_radial_uniform = createMeshCylindrical1D(Nr, Lr) 
# 1D radial nonuniform mesh
m_1d_radial_nonuniform = createMeshCylindrical1D(r) 
# 2D cartesian uniform mesh
m_2d_cartesian_uniform = createMesh2D(Nx, Ny, Lx, Ly) 
# 2D cartesian nonuniform mesh
m_2d_cartesian_nonuniform = createMesh2D(x, y) 
# 2D cylindrical uniform mesh
m_2d_cylindrical_uniform = createMeshCylindrical2D(Nr, Nz, Lr, Lz) 
# 2D cylindrical nonuniform mesh
m_2d_cylindrical_nonuniform = createMeshCylindrical2D(r, z) 
# 2D radial uniform mesh
m_2d_radial_uniform = createMeshRadial2D(Nr, Ntheta, Lr, Ltheta) 
# 2D radial nonuniform mesh
m_2d_radial_nonuniform = createMeshRadial2D(r, theta) 
# 3D cartesian uniform mesh
m_3d_cartesian_uniform = createMesh3D(Nx, Ny, Nz, Lx, Ly, Lz) 
# 3D cartesian nonuniform mesh
m_3d_cartesian_nonuniform = createMesh3D(x, y, z) 
# 3D cylindrical uniform mesh
m_3d_cylindrical_uniform = createMeshCylindrical3D(Nr, Ntheta, Nz, Lr, Ltheta, Lz) 
# 3D cylindrical nonuniform mesh
m_3d_cylindrical_nonuniform = createMeshCylindrical3D(r, theta, z)
println("Various mesh types in JFVM!")

Nx = 10
Lx = 1.0
m = createMesh1D(Nx, Lx)
BC = createBC(m)
BC.left.a[:]=BC.right.a[:]=0.0
BC.left.b[:]=BC.right.b[:]=1.0
BC.left.c[:]=1.0
BC.right.c[:]=0.0

c_init = 0.0 # initial value of the variable
c_old = createCellVariable(m, 0.0, BC)
c_old.value

D_val = 1.0 # value of the diffusion coefficient
D_cell = createCellVariable(m, D_val) # assigned to cells

# Harmonic average
D_face = harmonicMean(D_cell)
# Geometric average
D_face = geometricMean(D_cell)
# Arithmetic average
D_face = arithmeticMean(D_cell)
# UpwindMean, tvdMean, and linearMean can also be used for other purposes, e.g., linearizing a PDE

N_steps = 20 # number of time steps
dt= sqrt(Lx^2/D_val)/N_steps # time step

M_diff = diffusionTerm(D_face) # matrix of coefficient for diffusion term
(M_bc, RHS_bc)=boundaryConditionTerm(BC) # matrix of coefficient and RHS for the BC
for i =1:5
    (M_t, RHS_t)=transientTerm(c_old, dt, 1.0)
    M=M_t-M_diff+M_bc # add all the [sparse] matrices of coefficient
    RHS=RHS_bc+RHS_t # add all the RHS's together
    c_old = solveLinearPDE(m, M, RHS) # solve the PDE for the first time step, and replace the old value
end

Nx = 10
Lx = 1.0
m = createMesh1D(Nx, Lx)
BC = createBC(m)
BC.left.a[:]=BC.right.a[:]=0.0
BC.left.b[:]=BC.right.b[:]=1.0
BC.left.c[:]=1.0
BC.right.c[:]=0.0
c_init = 0.0 # initial value of the variable
c_old = createCellVariable(m, 0.0, BC)
D_val = 1.0 # value of the diffusion coefficient
D_cell = createCellVariable(m, D_val) # assigned to cells
# Harmonic average
D_face = harmonicMean(D_cell)
N_steps = 20 # number of time steps
dt= sqrt(Lx^2/D_val)/N_steps # time step
M_diff = diffusionTerm(D_face) # matrix of coefficient for diffusion term
(M_bc, RHS_bc)=boundaryConditionTerm(BC) # matrix of coefficient and RHS for the BC
for i =1:5
    (M_t, RHS_t)=transientTerm(c_old, dt, 1.0)
    M=M_t-M_diff+M_bc # add all the [sparse] matrices of coefficient
    RHS=RHS_bc+RHS_t # add all the RHS's together
    c_old = solveLinearPDE(m, M, RHS) # solve the PDE
    visualizeCells(c_old)
end

function trans_diff_dirichlet(m::MeshStructure)
    # a transient diffusion 
    BC = createBC(m)
    BC.left.a[:]=BC.right.a[:]=0.0
    BC.left.b[:]=BC.right.b[:]=1.0
    BC.left.c[:]=1.0
    BC.right.c[:]=0.0
    c_init = 0.0 # initial value of the variable
    c_old = createCellVariable(m, 0.0, BC)
    D_val = 1.0 # value of the diffusion coefficient
    D_cell = createCellVariable(m, D_val) # assigned to cells
    # Harmonic average
    D_face = harmonicMean(D_cell)
    N_steps = 20 # number of time steps
    dt= sqrt(m.facecenters.x[end]^2/D_val)/N_steps # time step
    M_diff = diffusionTerm(D_face) # matrix of coefficient for diffusion term
    (M_bc, RHS_bc)=boundaryConditionTerm(BC) # matrix of coefficient and RHS for the BC
    for i =1:5
        (M_t, RHS_t)=transientTerm(c_old, dt, 1.0)
        M=M_t-M_diff+M_bc # add all the [sparse] matrices of coefficient
        RHS=RHS_bc+RHS_t # add all the RHS's together
        c_old = solveLinearPDE(m, M, RHS) # solve the PDE
    end
    return c_old
end

Nx = 10
Ny = 15
Lx = 1.0
Ly = 2.0
m = createMesh2D(Nx, Ny, Lx, Ly)
c = trans_diff_dirichlet(m)
visualizeCells(c)

using PyPlot
Nx = 10
Ny = 15
Nz = 20
Lx = 1.0
Ly = 2.0
Lz = 3.0
m = createMesh3D(Nx, Ny, Nz, Lx, Ly, Lz)
c = trans_diff_dirichlet(m)
img = visualizeCells(c)
imshow(img)

Nr=10
Ntheta = 15
Lr = 1.0
Ltheta = 2*pi
m = createMeshRadial2D(Nr, Ntheta, Lr, Ltheta)
m.facecenters.x+=0.1
m.cellcenters.x+=0.1
c = trans_diff_dirichlet(m)
visualizeCells(c)

Nr=10
Ntheta = 15
Nz = 12
Lr = 1.0
Ltheta = 2*pi
Lz = 3.0
m = createMeshCylindrical3D(Nr, Ntheta, Nz, Lr, Ltheta, Lz)
m.facecenters.x+=0.1
m.cellcenters.x+=0.1
c = trans_diff_dirichlet(m)
img=visualizeCells(c)
imshow(img)