import numpy as np
import math
from create_basic_mesh import gen_block_mesh,get_tensor_product_conn
from sfepy.discrete.fem.meshio import UserMeshIO

def create_mesh_file():
    mesh = gen_block_mesh(np.array((0.001, 0.001)),
                          np.array((100,100)),
                          np.array((1.0, 1.0)),
                          name='')
    mesh.write('0.vtk', io = 'auto')

create_mesh_file()
#filename_mesh = UserMeshIO(mesh_hook)
#filename_mesh = "block.mesh"

filename_mesh = "0.vtk"

q = 1.6e-19
esi = 11.8*8.85e-12
ni2 = 1e20
v0 = 0.0259
dim = 0,001
xp = 0.0005
xn = 0.0005
na = 1e16
nd = 1e16

def doping_p(ts, coors, mode=None, **kwargs):
    x = coors[:,0]
    Na = np.where(x<0.0,na,0)
    return Na

def doping_n(ts, coors, mode=None, **kwargs):
    x = coors[:,0]
    Nd = np.where(x>0.0,nd,0) 
    return Nd

def charge_density(ts, coors, problem, equations=None, mode=None, **kwargs):
    if mode == 'qp':
        Na = doping_p(ts, coors, mode, **kwargs)
        Nd = doping_n(ts, coors, mode, **kwargs)
        V_values = problem.evaluate('ev_volume_integrate.i1.giunzione(V)',
                                    mode='qp', verbose=True)
        Na.shape = V_values.shape
        Nd.shape = V_values.shape
        ro=(-q/esi)*(Na*(-1+np.exp(-V_values/v0))+Nd*(1+np.exp(-V_values/v0)))
        #print V_values
        #print ro.shape, V_values.shape, Na.shape, Nd.shape
        ro.shape = (ro.shape[0]*ro.shape[1],1,1)
        return {'val': ro}

regions = {
    'giunzione' : 'all',
    'p' : ('vertices in (x <= 0.00000001 )', 'vertex'),
    'n' : ('vertices in (x >= 0.00099999 )', 'vertex'),
}

fields = {
    'potential': ('real',1,'giunzione',1),
}

variables = {
    'V' : ('unknown field', 'potential', 0),
    's' : ('test field',    'potential', 'V'),
   }

materials = {
    'coef': ({'val' : 1.0},),
    'coef1': 'charge_density',
    'Na': 'doping_p',
    'Nd': 'doping_n',
}

functions = {
    'charge_density' : (charge_density,),
    'doping_p' : (doping_p,),
    'doping_n' : (doping_n,),
}

ebcs = {
    'V1' : ('p',{'V.0': 0.0}),
}

integrals = {
    'i1' : 1
}

equations = {
    'potential' : """dw_laplace.i1.giunzione(coef.val, s, V) - 
    dw_volume_dot.i1.giunzione(coef1.val, s, V)
    =0""",
}


solvers = {
    'ls' : ('ls.scipy_direct', {}),
    'newton' : ('nls.newton', {
        'i_max'      : 100,
        'eps_a'      : 1e-10,
    }),
}

options = {
    'nls' : 'newton',
    'ls' : 'ls',
}