"""
Notes
-----

Important attributes of continuous (order > 0) :class:`Field` and
:class:`SurfaceField` instances:

- `vertex_remap` : `econn[:, :n_vertex] = vertex_remap[conn]`
- `vertex_remap_i` : `conn = vertex_remap_i[econn[:, :n_vertex]]`

where `conn` is the mesh vertex connectivity, `econn` is the
region-local field connectivity.
"""
from __future__ import absolute_import
import numpy as nm

from sfepy.base.base import output, get_default, assert_
from sfepy.base.base import Struct
from sfepy.discrete.common.fields import parse_shape, Field
from sfepy.discrete.fem.mesh import Mesh
from sfepy.discrete.fem.meshio import convert_complex_output
from sfepy.discrete.fem.utils import (extend_cell_data, prepare_remap,
                                      invert_remap, get_min_value)
from sfepy.discrete.fem.mappings import VolumeMapping, SurfaceMapping
from sfepy.discrete.fem.poly_spaces import PolySpace
from sfepy.discrete.fem.fe_surface import FESurface
from sfepy.discrete.integrals import Integral
from sfepy.discrete.fem.linearizer import (get_eval_dofs, get_eval_coors,
                                           create_output)
import six

def set_mesh_coors(domain, fields, coors, update_fields=False, actual=False,
                   clear_all=True, extra_dofs=False):
    if actual:
        if not hasattr(domain.mesh, 'coors_act'):
            domain.mesh.coors_act = nm.zeros_like(domain.mesh.coors)
        domain.mesh.coors_act[:] = coors[:domain.mesh.n_nod]
    else:
        domain.cmesh.coors[:] = coors[:domain.mesh.n_nod]

    if update_fields:
        for field in six.itervalues(fields):
            field.set_coors(coors, extra_dofs=extra_dofs)
            field.clear_mappings(clear_all=clear_all)

def eval_nodal_coors(coors, mesh_coors, region, poly_space, geom_poly_space,
                     econn, only_extra=True):
    """
    Compute coordinates of nodes corresponding to `poly_space`, given
    mesh coordinates and `geom_poly_space`.
    """
    if only_extra:
        iex = (poly_space.nts[:,0] > 0).nonzero()[0]
        if iex.shape[0] == 0: return

        qp_coors = poly_space.node_coors[iex, :]
        econn = econn[:, iex].copy()

    else:
        qp_coors = poly_space.node_coors

    ##
    # Evaluate geometry interpolation base functions in (extra) nodes.
    bf = geom_poly_space.eval_base(qp_coors)
    bf = bf[:,0,:].copy()

    ##
    # Evaluate extra coordinates with 'bf'.
    cmesh = region.domain.cmesh
    conn = cmesh.get_incident(0, region.cells, region.tdim)
    conn.shape = (econn.shape[0], -1)

    ecoors = nm.dot(bf, mesh_coors[conn])
    coors[econn] = nm.swapaxes(ecoors, 0, 1)

def _interp_to_faces(vertex_vals, bfs, faces):
    dim = vertex_vals.shape[1]
    n_face = faces.shape[0]
    n_qp = bfs.shape[0]

    faces_vals = nm.zeros((n_face, n_qp, dim), nm.float64)
    for ii, face in enumerate(faces):
        vals = vertex_vals[face,:dim]
        faces_vals[ii,:,:] = nm.dot(bfs[:,0,:], vals)

    return(faces_vals)

def get_eval_expression(expression,
                        fields, materials, variables,
                        functions=None, mode='eval', term_mode=None,
                        extra_args=None, verbose=True, kwargs=None):
    """
    Get the function for evaluating an expression given a list of elements,
    and reference element coordinates.
    """
    from sfepy.discrete.evaluate import eval_in_els_and_qp

    def _eval(iels, coors):
        val = eval_in_els_and_qp(expression, iels, coors,
                                 fields, materials, variables,
                                 functions=functions, mode=mode,
                                 term_mode=term_mode,
                                 extra_args=extra_args, verbose=verbose,
                                 kwargs=kwargs)
        return val[..., 0]

    return _eval

def create_expression_output(expression, name, primary_field_name,
                             fields, materials, variables,
                             functions=None, mode='eval', term_mode=None,
                             extra_args=None, verbose=True, kwargs=None,
                             min_level=0, max_level=1, eps=1e-4):
    """
    Create output mesh and data for the expression using the adaptive
    linearizer.

    Parameters
    ----------
    expression : str
        The expression to evaluate.
    name : str
        The name of the data.
    primary_field_name : str
        The name of field that defines the element groups and polynomial
        spaces.
    fields : dict
        The dictionary of fields used in `variables`.
    materials : Materials instance
        The materials used in the expression.
    variables : Variables instance
        The variables used in the expression.
    functions : Functions instance, optional
        The user functions for materials etc.
    mode : one of 'eval', 'el_avg', 'qp'
        The evaluation mode - 'qp' requests the values in quadrature points,
        'el_avg' element averages and 'eval' means integration over
        each term region.
    term_mode : str
        The term call mode - some terms support different call modes
        and depending on the call mode different values are
        returned.
    extra_args : dict, optional
        Extra arguments to be passed to terms in the expression.
    verbose : bool
        If False, reduce verbosity.
    kwargs : dict, optional
        The variables (dictionary of (variable name) : (Variable
        instance)) to be used in the expression.
    min_level : int
        The minimum required level of mesh refinement.
    max_level : int
        The maximum level of mesh refinement.
    eps : float
        The relative tolerance parameter of mesh adaptivity.

    Returns
    -------
    out : dict
        The output dictionary.
    """
    field = fields[primary_field_name]
    vertex_coors = field.coors[:field.n_vertex_dof, :]

    ps = field.poly_space
    gps = field.gel.poly_space
    vertex_conn = field.econn[:, :field.gel.n_vertex]

    eval_dofs = get_eval_expression(expression,
                                    fields, materials, variables,
                                    functions=functions,
                                    mode=mode, extra_args=extra_args,
                                    verbose=verbose, kwargs=kwargs)
    eval_coors = get_eval_coors(vertex_coors, vertex_conn, gps)

    (level, coors, conn,
     vdofs, mat_ids) = create_output(eval_dofs, eval_coors,
                                     vertex_conn.shape[0], ps,
                                     min_level=min_level,
                                     max_level=max_level, eps=eps)

    mesh = Mesh.from_data('linearized_mesh', coors, None, [conn], [mat_ids],
                          field.domain.mesh.descs)

    out = {}
    out[name] = Struct(name='output_data', mode='vertex',
                       data=vdofs, var_name=name, dofs=None,
                       mesh=mesh, level=level)

    out = convert_complex_output(out)

    return out

class FEField(Field):
    """
    Base class for finite element fields.

    Notes
    -----
    - interps and hence node_descs are per region (must have single
      geometry!)

    Field shape information:

    - ``shape`` - the shape of the base functions in a point
    - ``n_components`` - the number of DOFs per FE node
    - ``val_shape`` - the shape of field value (the product of DOFs and
      base functions) in a point
    """

    def __init__(self, name, dtype, shape, region, approx_order=1):
        """
        Create a finite element field.

        Parameters
        ----------
        name : str
            The field name.
        dtype : numpy.dtype
            The field data type: float64 or complex128.
        shape : int/tuple/str
            The field shape: 1 or (1,) or 'scalar', space dimension (2, or (2,)
            or 3 or (3,)) or 'vector', or a tuple. The field shape determines
            the shape of the FE base functions and is related to the number of
            components of variables and to the DOF per node count, depending
            on the field kind.
        region : Region
            The region where the field is defined.
        approx_order : int or tuple
            The FE approximation order. The tuple form is (order, has_bubble),
            e.g. (1, True) means order 1 with a bubble function.

        Notes
        -----
        Assumes one cell type for the whole region!
        """
        shape = 
parse_shape(shape, region.domain.shape.dim)