"""
Friction-slip model formulated as the implicit complementarity problem.

To integrate over a (dual) mesh, one needs:

* coordinates of element vertices
* element connectivity
* local base for each element
  * constant in each sub-triangle of the dual mesh

Data for each dual element:

* connectivity of its sub-triangles
* base directions t_1, t_2

Normal stresses:

* Assemble the rezidual and apply the LCBC operator described below.

Solution in \hat{V}_h^c:

* construct a restriction operator via LCBC just like in the no-penetration case
* use the substitution:
  u_1 = n_1 * w
  u_2 = n_2 * w
  u_3 = n_3 * w
  The new DOF is `w`.
* for the record, no-penetration does:
  w_1 = - (1 / n_1) * (u_2 * n_2 + u_3 * n_3)
  w_2 = u_2
  w_3 = u_3
"""
from sfepy.base.base import *
from sfepy.base.compat import unique
import sfepy.linalg as la
from sfepy.fem import Mesh, Domain, Field, Variables
from sfepy.fem.mappings import VolumeMapping, SurfaceMapping
from sfepy.fem.fe_surface import FESurface
from sfepy.fem.utils import compute_nodal_normals

def edge_data_to_output(coors, conn, e_sort, data):
    out = nm.zeros_like(coors)
    out[conn[e_sort,0]] = data
    return Struct(name='output_data',
                  mode='vertex', data=out,
                  dofs=None)

class DualMesh(Struct):
    """Dual mesh corresponding to a (surface) region."""

    def __init__(self, region):
        """
        Assume a single GeometryElement type in all groups, linear
        approximation.

        Works for one group only for the moment.
        """
        domain = region.domain

        self.dim = domain.shape.dim

        self.region = copy(region)
        self.region.setup_face_indices()

        self.mesh_coors = domain.mesh.coors

        # add_to_regions=True due to Field implementation shortcomings.
        omega = domain.create_region('Omega', 'all', add_to_regions=True)
        self.field = Field('displacements', nm.float64, (3,), omega, 1)

        self.gel = domain.geom_els.values()[0]
        self.sgel = self.gel.surface_facet

        face_key = 's%d' % self.sgel.n_vertex

        # Coordinate interpolation to face centres.
        self.ps = self.gel.interp.poly_spaces[face_key]
        centre = self.ps.node_coors.sum(axis=0) / self.ps.n_nod
        self.bf = self.ps.eval_base(centre[None,:])

        self.surfaces = surfaces = {}
        self.dual_surfaces = dual_surfaces = {}

        for ig, conn in enumerate(domain.mesh.conns):
            surface = FESurface(None, self.region, self.gel.faces, conn, ig)
            surfaces[ig] = surface

            dual_surface = self.describe_dual_surface(surface)
            dual_surfaces[ig] = dual_surface

    def describe_dual_surface(self, surface):
        n_fa, n_edge = surface.n_fa, self.sgel.n_edge

        mesh_coors = self.mesh_coors

        # Face centres.
        fcoors = mesh_coors[surface.econn]
        centre_coors = nm.dot(self.bf.squeeze(), fcoors)

        surface_coors = mesh_coors[surface.nodes]

        dual_coors = nm.r_[surface_coors, centre_coors]
        coor_offset = surface.nodes.shape[0]

        # Normals in primary mesh nodes.
        nodal_normals = compute_nodal_normals(surface.nodes, self.region,
                                              self.field)

        ee = surface.leconn[:,self.sgel.edges].copy()
        edges_per_face = ee.copy()
        sh = edges_per_face.shape
        ee.shape = edges_per_face.shape = (sh[0] * sh[1], sh[2])
        edges_per_face.sort(axis=1)

        eo = nm.empty((sh[0] * sh[1],), dtype=nm.object)
        eo[:] = [tuple(ii) for ii in edges_per_face]

        ueo, e_sort, e_id = 
unique(eo, return_index=True, return_inverse=True)