# Copyright 2016 Devsim LLC
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

from devsim import *
import sys
import math


def print_header(fh):
    '''
        Write header for backround mesh view
    '''
    fh.write('View "background mesh" {\n')

def print_footer(fh):
    '''
        Write footer for backround mesh view
    '''
    fh.write('};\n')

def get_edge_index(device, region):
    '''
    maps element edges to regular edges
    '''
    # now iterate over the edges of the element
    if "eindex" not in get_element_model_list(device=device, region=region):
        element_model(device=device, region=region,
                      name="eindex", equation="edge_index")
    eindex = get_element_model_values(
        device=device, region=region, name='eindex')
    eindex = [int(x) for x in eindex]
    return eindex

def get_node_index(device, region):
    '''
    maps head and tail nodes of from their edge index
    '''
    # identify all edges that need to be bisected
    # ultimately translated to an element
    if "node_index@n0" not in get_edge_model_list(device=device, region=region):
        edge_from_node_model(node_model="node_index", device=device, region=region)
    nindex = list(
        zip(
            [int(x) for x in get_edge_model_values(device=device, region=region, name="node_index@n0")],
            [int(x) for x in get_edge_model_values(device=device, region=region, name="node_index@n1")],
        )
    )
    return nindex

def calculate_clengths(device, region, model_values):
    '''
    calculate the characteristic lengths for each edge by bisecting the edge length
    '''
    clengths = list(get_edge_model_values(device=device, region=region, name="EdgeLength"))
    bisection_count = 0
    for i, v in enumerate(model_values):
        if v != 0:
            clengths[i] *= 0.5
            bisection_count += 1
    print("Edge Bisections: %d" % bisection_count)
    return clengths


def get_output_elements3(device, nindex, eindex, clengths, number_nodes, mincl, maxcl):
    '''
    gets the node indexes and the characterisic lengths for each element
    device : device we are operating on
    nindex : from get_node_index
    eindex : from get_edge_index
    clengths : from calculate_clengths
    number_nodes : number of nodes
    mincl : minimum characteristic length
    maxcl : maximum characteristic length
    '''
    # set upper limit to maxcl
    node_map = [maxcl] * number_nodes
    # get node indexes for each edge
    for i, n in enumerate(nindex):
        # clip minimum value to mincl
        v = max(clengths[i], mincl)
        for ni in n:
            node_map[ni] = min(node_map[ni], v) 

    dim = get_dimension(device=device)
    if dim == 2:
        # 3 edges per triangle
        skip = 3
    elif dim == 3:
        # 6 edges per tetrahedron
        skip = 6
    else:
        raise RuntimeError("Unhandled dimension %d" % dim)

    #break into a per element basis
    outputelements = []
    for i in range(0, len(eindex), skip):
        ndict = {}
        # mapping of element edge into an edge index
        for j in eindex[i:i+skip]:
            # mapping of edge index into a node index
            for k in nindex[j]:
                if k not in ndict:
                    ndict[k] = node_map[k]
        outputelements.append(tuple(ndict.items()))
    return outputelements

def print_elements(fh, device, region, elements):
    '''
    print background mesh triangles
    '''
    x = get_node_model_values(device=device, region=region, name="x")
    y = get_node_model_values(device=device, region=region, name="y")
    z = get_node_model_values(device=device, region=region, name="z")

    dim = get_dimension(device=device)
    if dim == 2:
        shape = "ST"
    elif dim == 3:
        shape = "SS"
    else:
        raise RuntimeError("Unhandled dimension %d" % dim)

    for e in elements:
        coords = []
        values = []
        for n, v in e:
            coords.extend((x[n], y[n], z[n]))
            values.append(v)
        coordstring = ", ".join([format(x, "1.15g") for x in coords])
        valuestring = ", ".join([format(x, "1.15g") for x in values])
        fh.write("%s(%s) {%s};\n" % (shape, coordstring, valuestring))

def refine_common(fh, device, region, model_values, mincl, maxcl):
    '''
    prints out the refined elements
    model_values : non-zero for edges to be bisected
    mincl : minimum characteristic length
    maxcl : maximum characteristic length
    '''
    clengths = calculate_clengths(device=device, region=region, model_values=model_values)

    eindex = get_edge_index(device, region)
    nindex = get_node_index(device, region)
    number_nodes = len(get_node_model_values(device=device, region=region, name="node_index"))


    outputelements = get_output_elements3(device=device, nindex=nindex, eindex=eindex, clengths=clengths, number_nodes=number_nodes, mincl=mincl, maxcl=maxcl)
    print_elements(fh=fh, device=device, region=region, elements=outputelements)

def get_oxide_model_values(device, region):
    '''
    returns a model for non-refinement
    '''
    test_model = [0.0] * len(get_edge_model_values(device=device, region=region, name="EdgeLength"))
    return test_model

def refine_oxide_region(fh, device, region, mincl, maxcl):
    '''
    refinement for oxide regions
    mincl : minimum characteristic length
    maxcl : maximum characteristic length
    '''
    # apply no refinement
    test_model = get_oxide_model_values(device=device, region=region)
    refine_common(fh=fh, device=device, region=region, model_values=test_model, mincl=mincl, maxcl=maxcl)

def get_silicon_model_values(device, region):
    '''
    returns a model for refinement of silicon regions
    '''
    # edge to node mapping (node0, node1)
    node_index = get_node_index(device=device, region=region)

    potential = get_node_model_values(device=device, region=region, name="Potential")
    test_model1 = [1 if abs(potential[x[0]]-potential[x[1]]) > 0.05 else 0 for x in node_index]

    electrons = get_node_model_values(device=device, region=region, name="Electrons")
    test_model2 = [1 if abs(math.log10(electrons[x[0]])-math.log10(electrons[x[1]])) > 1 else 0 for x in node_index]

    test_model = max_merge_lists((test_model1, test_model2))

    return test_model

def max_merge_lists(list_of_lists):
    '''
    gets the max of list of lists
    '''
    test_model = [max(x) for x in zip(*list_of_lists)]
    return test_model

def refine_silicon_region(fh, device, region, mincl, maxcl):
    '''
    refinement for silicon regions
    mincl : minimum characteristic length
    maxcl : maximum characteristic length
    '''

    test_model = get_silicon_model_values(device=device, region=region)

    refine_common(fh=fh, device=device, region=region, model_values=test_model, mincl=mincl, maxcl=maxcl)



