from __future__ import absolute_import
import numpy as nm

import sfepy.linalg as la
from sfepy.discrete.integrals import Integral
from sfepy.discrete import PolySpace
from six.moves import range

def prepare_remap(indices, n_full):
    """
    Prepare vector for remapping range `[0, n_full]` to its subset given
    by `indices`.
    """
    remap = nm.empty((n_full,), dtype=nm.int32)
    remap.fill(-1)
    remap[indices] = nm.arange(indices.shape[0], dtype=nm.int32)

    return remap

def invert_remap(remap):
    """
    Return the inverse of `remap`, i.e. a mapping from a sub-range
    indices to a full range, see :func:`prepare_remap()`.
    """
    if remap is not None:
        inverse = nm.where(remap >= 0)[0].astype(nm.int32)

    else:
        inverse = None

    return inverse

def prepare_translate(old_indices, new_indices):
    """
    Prepare vector for translating `old_indices` to `new_indices`.

    Returns
    -------
    translate : array
        The translation vector. Then `new_ar = translate[old_ar]`.
    """
    old_indices = nm.asarray(old_indices)
    new_indices = nm.asarray(new_indices)

    translate = nm.zeros(old_indices.max() + 1, dtype=new_indices.dtype)
    translate[old_indices] = new_indices

    return translate

def compute_nodal_normals(nodes, region, field, return_imap=False):
    """
    Nodal normals are computed by simple averaging of element normals of
    elements every node is contained in.
    """
    dim = region.dim

    field.domain.create_surface_group(region)
    field.setup_surface_data(region)

    # Custom integral with quadrature points in nodes.
    ps = PolySpace.any_from_args('', field.gel.surface_facet,
                                 field.approx_order)
    qp_coors = ps.node_coors
    # Unit normals -> weights = ones.
    qp_weights = nm.ones(qp_coors.shape[0], dtype=nm.float64)

    integral = 
Integral('aux', coors=qp_coors, weights=qp_weights)