# Source code for connectivity

    #
    ##
    ##  This file is part of pyFormex 1.0.5  (Sat Feb 16 10:40:32 CET 2019)
    ##  pyFormex is a tool for generating, manipulating and transforming 3D
    ##  geometrical models by sequences of mathematical operations.
    ##  Home page: http://pyformex.org
    ##  Project page:  http://savannah.nongnu.org/projects/pyformex/
    ##  Copyright 2004-2018 (C) Benedict Verhegghe (benedict.verhegghe@ugent.be)
    ##  Distributed under the GNU General Public License version 3 or later.
    ##
    ##  This program is free software: you can redistribute it and/or modify
    ##  it under the terms of the GNU General Public License as published by
    ##  the Free Software Foundation, either version 3 of the License, or
    ##  (at your option) any later version.
    ##
    ##  This program is distributed in the hope that it will be useful,
    ##  but WITHOUT ANY WARRANTY; without even the implied warranty of
    ##  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    ##  GNU General Public License for more details.
    ##
    ##  You should have received a copy of the GNU General Public License
    ##  along with this program.  If not, see http://www.gnu.org/licenses/.
    ##

    """A class and functions for handling nodal connectivity.

    This module defines a specialized array class for representing nodal
    connectivity. This is e.g. used in mesh models, where geometry is
    represented by a set of numbered points (nodes) and the geometric elements
    are described by refering to the node numbers.
    In a mesh model, points common to adjacent elements are unique, and
    adjacency of elements can easily be detected from common node numbers.
    """
    from __future__ import absolute_import, division, print_function

    from itertools import combinations
    import numpy as np

    from pyformex import arraytools as at
    from pyformex.adjacency import Adjacency

    ############################################################################
    ##
    ##   class Connectivity
    ##
    #########################
    #

    class Connectivity(np.ndarray):
        #
        # :DEV
        # Because we have a __new__ constructor and no __init__,
        # we have to put the signature of the object creation explicitely
        # in the first line of the docstring.
        #
        """Connectivity(data=[],dtyp=None,copy=False,nplex=0,eltype=None)

        A class for handling element to node connectivity.

        A connectivity object is a 2-dimensional integer array with all
        non-negative values. Each row of the array defines an element by listing
        the numbers of its lower entity types.
        A typical use is a :class:`~mesh.Mesh` object, where each element
        is defined in function of its nodes.
        While in a Mesh the word 'node' will normally refer to a geometrical
        point, here we will use 'node' for the lower entity whatever its nature
        is. It doesn't even have to be a geometrical entity.

        Note
        ----
        The current implementation limits a Connectivity object to numbers that
        are smaller than 2**31. That is however largely sufficient for all
        practical cases.

        In a row (element), the same node number may occur more than once, though
        usually all numbers in a row are different. Rows containing duplicate
        numbers are called `degenerate` elements.
        Rows containing the same node sets, albeit different permutations thereof,
        are called duplicates.

        Parameters
        ----------
        data: int array-like
            Data to initialize the Connectivity. The data should be 2-dim with
            shape ``(nelems,nplex)``, where ``nelems`` is the number of elements and
            ``nplex`` is the plexitude of the elements.
        dtyp: float datatype, optional
            It not provided, the datatype of ``data`` is used.
        copy: bool, optional
            If True, the data are copied. The default setting will try to use
            the original data if possible, e.g. if ``data`` is a correctly shaped
            and typed :class:`numpy.ndarray`.
        nplex: int, optional
            The plexitude of the data. This can be specified to force a check on
            the plexitude of the data, or to set the plexitude for an empty
            Connectivity.
            If an ``eltype`` is specified, the plexitude of the element type
            will override this value.
        eltype: str or :class:`elements.ElementType` subclass, optional
            The element type associated with the Connectivity. It can be either
            a subclass of:class:`elements.ElementType` or the ``name`` of such
            a subclass.
            If not provided, a non-typed Connectivity will result. If that is
            used to create a :class:`Mesh`, the proper element type will have to
            be specified at Mesh creation time. If the Connectivity will be used
            for other purposes, the element type may not be needed
            or not be important.

        Raises
        ------
        ValueError
            If ``nplex`` is provided and the specified ``data`` do not match the
            specified plexitude.

        Notes
        -----
        Empty Connectivities with ``nelems==0`` and ``nplex > 0`` can be useful,
        but a Connectivity with ``nplex==0`` generally is not.

        Examples
        --------
        >>> Connectivity([[0,1,2],[0,1,3],[0,3,2],[0,5,3]])
        Connectivity([[0, 1, 2],
                      [0, 1, 3],
                      [0, 3, 2],
                      [0, 5, 3]])

        >>> Connectivity(np.array([],dtype=at.Int).reshape(0,3))
        Connectivity([], shape=(0, 3))

        """

        _exclude_members_ = [ 'reorderNodes' ]

        def __new__(self,data=[],dtyp=None,copy=False,nplex=0,eltype=None):
            """Create a new Connectivity object."""
            if isinstance(data, Connectivity):
                if nplex == 0:
                    nplex = data.nplex()
                if eltype is None:
                    eltype = data.eltype

            if eltype is None:
                try:
                    eltype = data.eltype
                except:
                    eltype = None

            # Turn the data into an array, and copy if requested
            ar = np.array(data, dtype=dtyp, copy=copy)
            if ar.ndim < 2:
                if nplex > 0:
                    ar = ar.reshape(-1, nplex)
                else:
                    ar = ar.reshape(-1, 1)

            elif ar.ndim > 2:
                raise ValueError("Expected 2-dim data")

            # Make sure dtype is an int type
            if ar.dtype.kind != 'i':
                ar = ar.astype(at.Int)

            # Check values
            if ar.size > 0:
                if ar.max() >= 2**31 or (ar.min() < 0):
                    raise ValueError("Negative or too large positive value in data")
                if nplex > 0 and ar.shape[1] != nplex:
                    raise ValueError("Expected data of plexitude %s" % nplex)
            else:
                if nplex > 0:
                    ar = ar.reshape(0, nplex)

            # Transform 'subarr' from an ndarray to our new subclass.
            ar = ar.view(self)

            ## # Other data
            ar.eltype = eltype  # ! this may be a string!!!!!!!!!!!
            ar.inv = None   # inverse index
            return ar

        def __array_finalize__(self, obj):
            # reset the attributes from passed original object
            # all extra attributes added in __new__ should be reset here
            self.eltype = getattr(obj, 'eltype', None)
            self.inv = getattr(obj, 'inv', None)

        def __reduce__(self):
            """Reduce the object to a pickled state"""
            # Get the pickled ndarray state (as a list, so we can change it)
            object_state = list(np.ndarray.__reduce__(self))
            # Define our own state with the extra attributes we added
            subclass_state = (self.eltype, None)
            # Store both in place of the original ndarray state
            object_state[2] = (object_state[2], subclass_state)
            return tuple(object_state)

        def __setstate__(self, state):
            """Restore from pickled state"""
            # In __reduce__, we replaced ndarray's state with a tuple
            # of itself and our own state
            try:
                nd_state, own_state = state
                np.ndarray.__setstate__(self, nd_state)
                self.eltype, self.inv = own_state
            except:
                try:
                    # try to read legacy pickle format, which did not save
                    # the element type (rev < 2360)
                    np.ndarray.__setstate__(self, state)
                    print("WARNING: Connectivity was restored without element type!")
                except:
                    print("I could not unpickle the Connectivity, neither in old not new format")
                    raise

        def __repr__(self):
            """String representation of a Coords

            Examples
            --------
            >>> Connectivity([[0,1,2],[0,1,3]],eltype='line3')
            Connectivity([[0, 1, 2],
                          [0, 1, 3]], eltype='line3')

            """
            res = np.ndarray.__repr__(self)
            if self.dtype == at.Int:
                res = res.replace(', dtype=int32','')
            if self.eltype is not None:
                res = res.replace(')',", eltype='%s')" % self.eltype)
            return res

        def nelems(self):
            """Return the number of elements in the Connectivity table.

            Returns
            -------
            int
                The number of rows in the table.

            Examples
            --------
            >>> Connectivity([[0,1,2],[0,1,3],[0,3,2],[0,5,3]]).nelems()
            4
            """
            return self.shape[0]

        def maxnodes(self):
            """Return an upper limit for number of nodes in the Connectivity.

            Returns
            -------
            int
                The highest node number plus one.

            See Also
            --------
            nnodes: the actual number of nodes in the table

            Examples
            --------
            >>> Connectivity([[0,1,2],[0,1,3],[0,3,2],[0,5,3]]).maxnodes()
            6

            """
            return int(self.max() + 1)

        def nnodes(self):
            """Return the actual number of nodes in the Connectivity.

            This returns the count of the unique node numbers.

            See Also
            --------
            maxnodes: the highest node number + 1

            Examples
            --------
            >>> Connectivity([[0,1,2],[0,1,3],[0,3,2],[0,5,3]]).nnodes()
            5
            """
            return np.unique(self).shape[0]

        def nplex(self):
            """Return the plexitude of the elements in the Connectivity table.

            Examples
            --------
            >>> Connectivity([[0,1,2],[0,1,3],[0,3,2],[0,5,3]]).nplex()
            3
            """
            return self.shape[1]

        def report(self):
            """Format a Connectivity table"""
            s = "Connectivity %s, eltype=%s" % (self.shape, self.eltype)
            return s + '\n' + np.ndarray.__str__(self)

    ############### Detecting degenerates and duplicates ##############

        def testDegenerate(self):
            """Flag the degenerate elements (rows).

            A degenerate element is a row which contains at least two
            equal values.

            Returns
            -------
            bool array
                 A 1-dim bool array with length ``self.nelems()``, holding
                 True values for the degenerate rows.

            Examples
            --------
            >>> Connectivity([[0,1,2],[0,1,1],[0,3,2]]).testDegenerate()
            array([False,  True, False], dtype=bool)

            """
            srt = np.asarray(self.copy())
            srt.sort(axis=1)
            return (srt[:, :-1] == srt[:, 1:]).any(axis=1)

        def listDegenerate(self):
            """Return a list with the numbers of the degenerate elements.

            Returns
            -------
            int array
                A 1-dim int array holding the row indices of the
                degenerate elements.

            Examples
            --------
            >>> Connectivity([[0,1,2],[0,1,1],[0,3,2]]).listDegenerate()
            array([1])

            """
            return np.arange(self.nelems())[self.testDegenerate()]

        def listNonDegenerate(self):
            """Return a list with the numbers of the non-degenerate elements.

            Returns
            -------
            int array
                A 1-dim int array holding the row indices of the
                non-degenerate elements.

            Examples
            --------
            >>> Connectivity([[0,1,2],[0,1,1],[0,3,2]]).listNonDegenerate()
            array([0, 2])

            """
            return np.arange(self.nelems())[~self.testDegenerate()]

        def removeDegenerate(self):
            """Remove the degenerate elements from a Connectivity table.

            Returns
            -------
            Connectivity
                A Connectivity object with the degenerate elements removed.

            Examples
            --------
            >>> Connectivity([[0,1,2],[0,1,1],[0,3,2]]).removeDegenerate()
            Connectivity([[0, 1, 2],
                          [0, 3, 2]])

            """
            return self[~self.testDegenerate()]

        def reduceDegenerate(self,target=None,return_indices=False):
            """Reduce degenerate elements to lower plexitude elements.

            This will try to reduce the degenerate elements of the Connectivity
            to a lower plexitude. This is only possible if an element type
            was set in the Connectivity. This function uses the data of
            the Element database in :mod:`elements` to reduce certain element
            types to other types depending on the kind of degeneracy.

            Parameters
            ----------
            target: str, optional
                Target element name. If provided, only reductions to that element
                type are performed. Else, all the target element types for which
                a reduction scheme is available, will be tried.
            return_indices: bool, optional
                If True, also returns the indices of the elements in the
                input Connectivity that are contained in each of the parts.

            Returns
            -------
            conn: list of Connectivity objects
                A list of Connectivities of which the first one contains
                the originally non-degenerate elements, the next one(s) contain
                the reduced elements (per reduced element type) and the last one
                contains elements that could not be reduced (may be absent).
                In the following cases a list with only the original is returned:

                - there are no degenerate elements in the Connectivity;
                - the Connectivity does not have an element type;
                - the element type of the Connectivity does not have any reduction
                  scheme defined.
            ind: list of indices, optional
                Only returned if ``return_indices`` is True: the indices of the
                elements of ``self`` contained in each item in ``conn``.

            Note
            ----
            If the Connectivity is part of a :class:`~mesh.Mesh`, you should
            use the :meth:`mesh.Mesh.splitDegenerate` method instead, as that
            will preserve the property numbers in the resulting Meshes.

            The resulting reduced Connectivities may still be degenerate for
            their own element type.

            See Also
            --------
            mesh.Mesh.splitDegenerate: split mesh in non-degenerate and reduced

            Examples
            --------
            >>> Connectivity([[0,1,2],[0,1,3]],eltype='line3').reduceDegenerate()
            [Connectivity([[0, 1, 2],
                          [0, 1, 3]], eltype='line3')]
            >>> C = Connectivity([[0,1,2],[0,1,1],[0,3,2]],eltype='line3')
            >>> C.reduceDegenerate()
            [Connectivity([[0, 1, 2],
                          [0, 3, 2]], eltype='line3'), \
            Connectivity([[0, 1]], eltype='line2')]
            >>> C = Connectivity([[0,1,2],[0,1,1],[0,3,2],[1,1,1],[0,0,2]],eltype='line3')
            >>> C.reduceDegenerate()
            [Connectivity([[0, 1, 2],
                          [0, 3, 2]], eltype='line3'), \
            Connectivity([[1, 1],
                          [0, 2],
                          [0, 1]], eltype='line2')]

            >>> C.reduceDegenerate(return_indices=True)
            ([Connectivity([[0, 1, 2],
                          [0, 3, 2]], eltype='line3'), Connectivity([[1, 1],
                          [0, 2],
                          [0, 1]], eltype='line2')], \
            [array([0, 2]), array([3, 4, 1])])
            """
            from pyformex.elements import elementType
            if self.eltype is None:
               return [ self ]

            eltype = elementType(self.eltype)
            if not hasattr(eltype, 'degenerate'):
                return [ self ]

            degen = self.testDegenerate()
            if not degen.any():
                return [ self ]

            # get all reductions for this eltype
            strategies = eltype.degenerate

            # if target, keep only those leading to target
            if target is not None:
                s = strategies.get(target, [])
                if s:
                    strategies = {target:s}
                else:
                    strategies = {}

            if not strategies:
                return [self]

            # non-degenerate: keep
            conn = [ self[~degen] ]
            ind = [ np.where(~degen)[0] ]

            # degenerate: reduce
            e = self[degen]
            i = np.where(degen)[0]

            for totype in strategies:

                ids = np.array([],dtype=at.Int)
                elems = []

                for conditions, selector in strategies[totype]:
                    cond = np.array(conditions)
                    w = (e[:, cond[:, 0]] == e[:, cond[:, 1]]).all(axis=1)
                    sel = np.where(w)[0]
                    if len(sel) > 0:
                        elems.append(e[sel][:, selector])
                        ids = np.concatenate([ids,i[sel]])
                        # remove the reduced elems from e
                        e = e[~w]
                        i = i[~w]
                        if e.nelems() == 0:
                            break

                if elems:
                    elems = np.concatenate(elems)
                    conn.append(Connectivity(elems, eltype=totype))
                    ind.append(ids)

                if e.nelems() == 0:
                    break

            if e.nelems() > 0:
                conn.append(e)
                ind.append(ids)

            if return_indices:
                return conn, ind
            else:
                return conn

        def findDuplicate(self,permutations='all'):
            """Find duplicate rows in the Connectivity.

            Parameters
            ----------
            permutations: str
                Defines which permutations of the row data are allowed while still
                considering the rows equal. Possible values are:

                - 'none': no permutations are allowed: rows must match the same date
                  at the same positions.
                - 'roll': rolling is allowed. Rows that can be transformed into
                  each other by rolling are considered equal;
                - 'all': any permutation of the same data will be considered an
                  equal row. This is the default.

            Returns
            -------
            V: :class:`~varray.Varray`
                A Varray where each row contains a list of the row numbers
                from a that are considered equal. The entries in each row are
                sorted and the rows are sorted according to their first element.

            Notes
            -----
            This is like :func:`arraytools.equalRows` but has a different
            default value for ``permutations``.

            Examples
            --------
            >>> C = Connectivity([[0,1,2],[0,1,3],[0,1,2],[2,0,1],[2,1,0]])
            >>> C.findDuplicate()
            Varray([[0, 2, 3, 4], [1]])
            >>> C.findDuplicate(permutations='roll')
            Varray([[0, 2, 3], [1], [4]])
            >>> C.findDuplicate(permutations='none')
            Varray([[0, 2], [1], [3], [4]])

            """
            return at.equalRows(self,permutations=permutations)

        def listDuplicate(self,permutations='all'):
            """Return a list with the numbers of the duplicate elements.

            Returns
            -------
            1-dim int array
                The indices of the unique rows in the Connectivity array.

            Examples
            --------
            >>> C = Connectivity([[0,1,2],[0,1,3],[0,1,2],[2,0,1],[2,1,0]])
            >>> C.listDuplicate()
            array([2, 3, 4])
            >>> C.listDuplicate(permutations='roll')
            array([2, 3])
            >>> C.listDuplicate(permutations='none')
            array([2])

            """
            ind, ok = at.findEqualRows(self,permutations=permutations)
            return np.sort(ind[~ok])

        def listUnique(self,permutations='all'):
            """Return a list with the numbers of the unique elements.

            Returns
            -------
            1-dim int array
                The indices of the unique rows in the Connectivity array.

            See Also
            --------
            findDuplicate: find duplicate rows
            listDuplicate: list duplicate rows
            removeDuplicate: remove duplicate rows

            Examples
            --------
            >>> C = Connectivity([[0,1,2],[0,1,3],[0,1,2],[2,0,1],[2,1,0]])
            >>> C.listUnique()
            array([0, 1])
            >>> C.listUnique(permutations='roll')
            array([0, 1, 4])
            >>> C.listUnique(permutations='none')
            array([0, 1, 3, 4])

            """
            return at.uniqueRows(self,permutations=permutations)

        def removeDuplicate(self,permutations='all'):
            """Remove duplicate elements from a Connectivity list.

            By default, duplicates are elements that consist of the same set of
            nodes, in any particular order. Setting permutations to 'none'
            will only remove the duplicate rows that have matching values at
            matching positions.

            Returns
            -------
            Connectivity
                A new Connectivity with the duplicate elements removed.

            Examples
            --------
            >>> C = Connectivity([[0,1,2],[0,1,3],[0,1,2],[2,0,1],[2,1,0]])
            >>> C.removeDuplicate()
            Connectivity([[0, 1, 2],
                          [0, 1, 3]])
            >>> C.removeDuplicate(permutations='roll')
            Connectivity([[0, 1, 2],
                          [0, 1, 3],
                          [2, 1, 0]])
            >>> C.removeDuplicate(permutations='none')
            Connectivity([[0, 1, 2],
                          [0, 1, 3],
                          [2, 0, 1],
                          [2, 1, 0]])

            """
            return self[self.listUnique(permutations)]

        def reorder(self,order='nodes'):
            """Reorder the elements of a Connectivity in a specified order.

            This does not actually reorder the elements itself, but returns
            an index with the order of the rows (elements) in the Connectivity
            table that meets the specified ordering requirements.

            Parameters
            ----------
            order: str or list of ints
                Specifies how to reorder the elements. It is either one
                of the special string values defined below, or else it is an index
                with length equal to the number of elements. The index should be
                a permutation of the numbers in ``range(self.nelems()``. Each value
                gives the number of the old element that should be placed at
                this position. Thus, the order values are the old element numbers
                on the position of the new element number.

                ``order`` can also take one of the following predefined values,
                resulting in the corresponding renumbering scheme being generated:

                - 'nodes': the elements are renumbered in order of their appearance
                  in the inverse index, i.e. first are the elements connected to
                  node 0, then the as yet unlisted elements connected to node 1, etc.
                - 'random': the elements are randomly renumbered.
                - 'reverse': the elements are renumbered in reverse order.

           	Returns
            -------
            1-dim int array
                Int array with a permutation of ``arange(self.nelems()``, such that
                taking the elements in this order will produce a Connectivity
                reordered as requested.
                In case an explicit order was specified as input, this order is
                returned after checking that it is indeed a permutation of
                ``range(self.nelems()``.

            Examples
            --------
            >>> A = Connectivity([[1,2],[2,3],[3,0],[0,1]])
            >>> A[A.reorder('reverse')]
            Connectivity([[0, 1],
                          [3, 0],
                          [2, 3],
                          [1, 2]])
            >>> A[A.reorder('nodes')]
            Connectivity([[0, 1],
                          [3, 0],
                          [1, 2],
                          [2, 3]])
            >>> A[A.reorder([2,3,0,1])]
            Connectivity([[3, 0],
                          [0, 1],
                          [1, 2],
                          [2, 3]])

            """
            if order == 'nodes':
                a = np.sort(self, axis=-1)  # first sort rows
                order = at.sortByColumns(a)
            elif order == 'reverse':
                order = np.arange(self.nelems()-1, -1, -1)
            elif order == 'random':
                order = random.permutation(self.nelems())
            else:
                order = np.asarray(order)
                if not (order.dtype.kind == 'i' and \
                        (np.sort(order) == np.arange(order.size)).all()):
                    raise ValueError("order should be a permutation of range(%s)" % self.nelems())
            return order

        def renumber(self,start=0):
            """Renumber the nodes to a consecutive integer range.

            The node numbers in the table are changed thus that they
            form a consecutive integer range starting from the specified
            value.

            Parameters
            ----------
            start: int
                Lowest node number to be used in the renumbered Connectivity.

            Returns
            -------
            elems: Connectivity
                The renumbered Connectivity
            oldnrs: 1-dim int array
                The sorted list of unique (old) node numbers. The new
                node numbers are assigned in order of increasing old node numbers,
                thus the old node number for new node number ``i`` can be found
                at position ``i - start``.

            Examples
            --------
            >>> e,n = Connectivity([[0,2],[1,4],[4,2]]).renumber(7)
            >>> print(e)
            [[ 7  9]
             [ 8 10]
             [10  9]]
            >>> print(n)
            [0 1 2 4]

            Find the old node number of new node 10
            >>> n[10-7]
            4
            """
            nodes = np.asarray(np.unique(self))
            if nodes.size == 0:
                elems = self

            else:
                old = np.arange(nodes.max()+1)
                if nodes.shape[0] == old.shape[0]:
                    # we have a consecutive range
                    if nodes[0] == start:
                        # numbering is ok, keep
                        elems = self
                    else:
                        # add the correct offset
                        elems = self + (start-nodes[0])
                else:
                    # need to renumber
                    elems = at.inverseUniqueIndex(nodes)[self] + start
                    elems = Connectivity(elems, eltype=self.eltype)

            return elems, nodes

        def inverse(self,expand=True):
            """Return the inverse index of a Connectivity table.

            Returns
            -------
            int array
                The inverse index of the Connectivity, as computed
                by :func:`arraytools.inverseIndex`.

            Examples
            --------
            >>> Connectivity([[0,1,2],[0,1,4],[0,4,2]]).inverse()
            array([[ 0,  1,  2],
                   [-1,  0,  1],
                   [-1,  0,  2],
                   [-1, -1, -1],
                   [-1,  1,  2]], dtype=int32)
            >>> Connectivity([[0,1,2],[0,1,4],[0,4,2]]).inverse(expand=False)
            Varray([[0, 1, 2], [0, 1], [0, 2], [], [1, 2]])
            """
            if self.inv is None or self.flags.writeable:
                if self.size > 0:
                    self.inv = at.inverseIndex(self,expand=expand)
                else:
                    self.inv = Connectivity()
            return self.inv

        def nParents(self):
            """Return the number of elements connected to each node.

            Returns
            -------
            1-dim int array
                The number of elements connected to each node. The length of
                the array is equal to the highest node number + 1.
                Unused node numbers will have a count of zero.

            Examples
            --------
            >>> Connectivity([[0,1,2],[0,1,4],[0,4,2]]).nParents()
            array([3, 2, 2, 0, 2])
            """
            r = self.inverse()
            return (r>=0).sum(axis=1)

        def connectedTo(self,nodes,return_ncon=False):
            """Check if the elements are connected to the specified nodes.

            Parameters
            ----------
            nodes: int or int array-like
                One or more node numbers to check for connections in the table.
            return_ncon: bool, optional
                If True, also return the number of connections for each element.

            Returns
            -------
            connections: int array
                The numbers of the elements that contain at least one of the
                specified nodes.
            ncon: int array, optional
                The number of connections for each connected element.
                This is only provided if ``return_ncon`` is True.

            Examples
            --------
            >>> A = Connectivity([[0,1,2],[0,1,3],[0,3,2],[1,2,3]])
            >>> print(A.connectedTo(2))
            [0 2 3]
            >>> A.connectedTo([0,1,3],True)
            (array([0, 1, 2, 3], dtype=int32), array([2, 3, 2, 2]))

            """
            nodes = at.checkArray1D(nodes, kind='i')
            nodes = np.intersect1d(nodes, self) #remove unconnected nodes
            inv = self.inverse()
            ad = inv[nodes]
            ad = ad[ad>=0]
            # We now have a list of all individual attachements to any of the nodes,
            # identified by the element number. We count them per element.
            m, u = at.multiplicity(ad)
            if return_ncon:
                return u, m
            else:
                return u

        def hits(self, nodes):
            """Count the nodes from a list connected to the elements.

            Parameters
            ----------
            nodes: int or list of ints
                One or more node numbers.

            Returns
            -------
            int array (nelems,)
                An int array holding the number of nodes from the specified
                input that are contained in each of the elements.

            Notes
            -----
            This information can also be got from meth:`connectedTo`.
            This method however expands the results to the full element set,
            making it apt for use in selector expressions like
            ``self[self.hits(nodes) >= 2]``.

            Examples
            --------
            >>> A = Connectivity([[0,1,2],[0,1,3],[0,3,2],[1,2,3]])
            >>> A.hits(2)
            array([1, 0, 1, 1])
            >>> A.hits([0,1,3])
            array([2, 3, 2, 2])
            """
            u, m = self.connectedTo(nodes, True)
            res = np.zeros(self.shape[0], dtype=m.dtype)
            res[u] = m
            return res

        def adjacency(self,kind='e',mask=None):
            """Create a table of adjacent items.

            This creates an element adjacency table or node adjacency table
            An element `i` is said to be adjacent to element `j`, if the two
            elements have at least one common node.
            A node `i` is said to be adjacent to node `j`, if there is at least
            one element containing both nodes.

            Parameters
            ----------
            kind: 'e' or 'n'
                Select element ('e') or node (n') adjacency table. Default is
                element adjacency.
            mask: bool array or int index, optional
                Node selector. If provided (with ``kind=='e'``) this defines
                by a bool flag array or int index numbers the list of nodes
                that are to be considered connectors between elements. The default
                is to consider all nodes as connectors.

                This option is only useful in the case `kind` == 'e'. If you want to
                use an element mask for the 'n' case, just apply the (element) mask
                beforehand by using ``self[mask].adjacency('n')``.

            Returns
            -------
            :class:`~adjacency.Adjacency` object
                An Adjacency array with shape (nr,nc), where row `i` holds
                a sorted list of all the items that are adjacent to item `i`,
                padded with -1 values to create an equal list length for all items.

            Examples
            --------
            >>> Connectivity([[0,1],[0,2],[1,3],[0,5]]).adjacency('e')
            Adjacency([[ 1,  2,  3],
                       [-1,  0,  3],
                       [-1, -1,  0],
                       [-1,  0,  1]], dtype=int32)
            >>> Connectivity([[0,1],[0,2],[1,3],[0,5]]).adjacency('e',mask=[1,2,3,5])
            Adjacency([[ 2],
                       [-1],
                       [ 0],
                       [-1]], dtype=int32)
            >>> Connectivity([[0,1],[0,2],[1,3],[0,5]]).adjacency('n')
            Adjacency([[ 1,  2,  5],
                       [-1,  0,  3],
                       [-1, -1,  0],
                       [-1, -1,  1],
                       [-1, -1, -1],
                       [-1, -1,  0]])
            >>> Connectivity([[0,1,2],[0,1,3],[2,4,5]]).adjacency('n')
            Adjacency([[-1,  1,  2,  3],
                       [-1,  0,  2,  3],
                       [ 0,  1,  4,  5],
                       [-1, -1,  0,  1],
                       [-1, -1,  2,  5],
                       [-1, -1,  2,  4]])
            >>> Connectivity([[0,1,2],[0,1,3],[2,4,5]])[[0,2]].adjacency('n')
            Adjacency([[-1, -1,  1,  2],
                       [-1, -1,  0,  2],
                       [ 0,  1,  4,  5],
                       [-1, -1, -1, -1],
                       [-1, -1,  2,  5],
                       [-1, -1,  2,  4]])
            """
            inv = self.inverse()
            if kind == 'e':
                if mask is not None:
                    mask = at.complement(mask, inv.shape[0])
                    inv[mask] = -1
                if self.size <= 4000000:
                    # do in one step
                    adj = _elem_adj(inv,self,True)
                else:
                    # use multiprocessing
                    # it saves a bit of memory, but not too much
                    # as we process in parallel
                    adj = _elem_adj_multi(inv,self,nproc=4)
                    maxcols = max([a.shape[1] for a in adj])
                    adj = [at.growAxis(a,maxcols-a.shape[1],axis=1,fill=-1)
                           for a in adj]
                    adj = np.concatenate(adj,axis=0)
                    adj = Adjacency(adj)

            elif kind == 'n':
                adj = np.concatenate([np.where(inv>=0, self[:, i][inv], inv)
                                      for i in range(self.nplex())], axis=1)
                adj = Adjacency(adj)
            else:
                raise ValueError("kind should be 'e' or 'n', got %s" % str(kind))
            return adj

        def adjacentElements(self,els,mask=None):
            """Compute adjacent elements.

            This creates an element adjacency table or node adjacency table
            An element `i` is said to be adjacent to element `j`, if the two
            elements have at least one common node.
            A node `i` is said to be adjacent to node `j`, if there is at least
            one element containing both nodes.

            Parameters
            ----------
            else: int or list of ints
                The element number(s) for which to compute the adjacent elements
            mask: bool array or int index, optional
                Node selector. If provided (with ``kind=='e'``) this defines
                by a bool flag array or int index numbers the list of nodes
                that are to be considered connectors between elements. The default
                is to consider all nodes as connectors.

                This option is only useful in the case `kind` == 'e'. If you want to
                use an element mask for the 'n' case, just apply the (element) mask
                beforehand by using ``self[mask].adjacency('n')``.

            Returns
            -------
            :class:`~adjacency.Adjacency` object
                An Adjacency array with shape (nr,nc), where row `i` holds
                a sorted list of all the items that are adjacent to item `i`,
                padded with -1 values to create an equal list length for all items.

            Examples
            --------
            >>> Connectivity([[0,1],[0,2],[1,3],[0,5]]).adjacentElements([0,1,2,3])
            array([[ 1,  2,  3],
                   [-1,  0,  3],
                   [-1, -1,  0],
                   [-1,  0,  1]], dtype=int32)
            >>> Connectivity([[0,1],[0,2],[1,3],[0,5]]).adjacentElements([0,1,2])
            array([[ 1,  2,  3],
                   [-1,  0,  3],
                   [-1, -1,  0]], dtype=int32)
            >>> Connectivity([[0,1],[0,2],[1,3],[0,5]]).adjacentElements([1,2,3])
            array([[ 0,  3],
                   [-1,  0],
                   [ 0,  1]], dtype=int32)
            >>> Connectivity([[0,1],[0,2],[1,3],[0,5]]).adjacentElements([0,2])
            array([[ 1,  2,  3],
                   [-1, -1,  0]], dtype=int32)
            >>> Connectivity([[0,1],[0,2],[1,3],[0,5]]).adjacentElements([2])
            array([[0]], dtype=int32)
            >>> Connectivity([[0,1],[0,2],[1,3],[0,5]]).adjacentElements(1)
            array([[0, 3]], dtype=int32)
            """
            els = at.checkArray1D(els,kind='i')
            nels = els.shape[0]
            if nels <= 0:
                return np.array([],dtype=at.Int)
            inv = self.inverse()
            #print(inv)
            if mask is not None:
                mask = at.complement(mask, inv.shape[0])
                inv[mask] = -1
            #print(inv[self[els]])
            adj = inv[self[els]].reshape((nels, -1))
            adj[adj == els.reshape(nels, -1)] = -1 # remove the element itself
            adj.sort(axis=1)
            adj[np.where(adj[:, :-1] == adj[:, 1:])] = -1 #remove duplicate items
            adj.sort(axis=1)
            maxc = adj.max(axis=0) # find maximum per column
            adj = adj[:, maxc>=0]   # retain columns with non-negative maximum
            return adj

        ### frontal methods ###

        # TODO: This is still much slower that Adjacency.frontwalk
        #   Maybe we should just remove this (or else implement in C)

        def frontGenerator(self,startat=0,frontinc=1,partinc=1):
            """Generator function returning the frontal elements.

            This is a generator function and is normally not used directly,
            but via the :meth:`frontWalk` method.

            Parameters: see :meth:`frontWalk`.

            Returns
            -------
            int array
                Int array with a value for each element. On the initial call,
                all values are -1, except for the elements in the initial front,
                which get a value 0. At each call a new front is created with
                all the elements that are connected to any of the
                current front and which have not yet been visited. The new front
                elements get a value equal to the last front's value plus the
                ``frontinc``. If the front becomes empty and a new starting front is
                created, the front value is extra incremented with ``partinc``.

            Examples
            --------
            >>> C = Connectivity([[2,8,7],[2,3,8],[3,9,8],[4,10,9],[5,6,11],
            ...                   [6,12,11]])
            >>> C.adjacentElements([0])
            array([[1, 2]], dtype=int32)
            >>> for p in C.frontGenerator(): print(p)
            [ 0 -1 -1 -1 -1 -1]
            [ 0  1  1 -1 -1 -1]
            [ 0  1  1  2 -1 -1]
            [ 0  1  1  2  4 -1]
            [0 1 1 2 4 5]
            >>> A = C.adjacency()
            >>> for p in A.frontGenerator(): print(p)
            [ 0 -1 -1 -1 -1 -1]
            [ 0  1  1 -1 -1 -1]
            [ 0  1  1  2 -1 -1]
            [ 0  1  1  2  4 -1]
            [0 1 1 2 4 5]
            """
            if self.nelems() <= 0:
                return

            # Initialize result
            p = -np.ones((self.nelems()), dtype=at.Int)

            # Compute inverse once, then lock array
            self.inv = self.inverse()
            self.flags.writeable = False

            # Remember current elements front
            elems = np.clip(np.asarray(startat), 0, self.nelems())
            prop = 0
            while elems.size > 0:
                # Store prop value for current elems
                #at.printar("Elems ",elems)
                p[elems] = prop
                yield p

                prop += frontinc

                # Determine adjacent elements
                #nodes = np.unique(np.asarray(self[elems]))
                #elems = self.connectedTo(nodes)
                elems = np.unique(self.adjacentElements(elems))
                while elems[0] < 0: # There should only be one -1
                    elems = elems[1:]
                # Remove already done
                elems = np.setdiff1d(elems, np.where(p>=0)[0])
                if elems.size > 0:
                    continue

                # No more elements in this part: start a new one
                elems = np.where(p<0)[0]
                if elems.size > 0:
                    # Start a new part
                    elems = elems[[0]]
                    prop += partinc

            # Free array
            self.inv = None
            self.flags.writeable = True

        def frontWalk(self,startat=0,frontinc=1,partinc=1,maxval=-1):
            """Walks through the elements by their node front.

            A frontal walk is executed starting from the given element(s).
            A number of steps is executed, each step advancing the front
            over a given number of single pass increments. The step number at
            which an element is reached is recorded and returned.

            Parameters
            ----------
            startat: int or list of ints
                Initial element number(s) in the front.
            frontinc: int
                Increment for the front number on each frontal step.
            partinc: int
                Increment for the front number when the front gets empty and
                a new part is started.
            maxval: int
                Maximum frontal value. If negative (default) the walk will
                continue until all elements have been reached. If non-negative,
                walking will stop as soon as the frontal value reaches this
                maximum.

            Returns
            -------
            int array
                An array of ints specifying for each element in which step
                the element was reached by the walker.

            Examples
            --------
            >>> C = Connectivity([[2,8,7],[2,3,8],[3,9,8],[4,10,9],[5,6,11],
            ...                   [6,12,11]])
            >>> print(C.frontWalk())
            [0 1 1 2 4 5]
            """
            for p in self.frontGenerator(startat=startat, frontinc=frontinc, partinc=partinc):
                if maxval >= 0:
                    if p.max() >= maxval:
                        break
            return p

        def front(self,startat=0,add=False):
            """Returns the elements of the first node front.

            Parameters
            ----------
            startat: int or list od ints
                Element number(s) or a list of element numbers. The list of
                elements to find the next front for.
            add: bool, optional
                If True, the `startat` elements wil be included in the
                return value. The default (False) will only return the elements
                in the next front line.

            Returns
            -------
            int array
                A list of the elements that are connected to any of the nodes
                that are part of the startat elements.

            Notes
            -----
            This is equivalent to the first step of a :func:`frontWalk`
            with the same startat elements, and could thus also be obtained from
            ``where(self.frontWalk(startat,maxval=1) == 1)[0]``.

            Here however another implementation is used, which is more efficient
            for very large models: it avoids the creation of the large array as
            returned by frontWalk.

            Examples
            --------
            >>> C = Connectivity([[2,8,7],[2,3,8],[3,9,8],[4,10,9],[5,6,11],
            ...                   [6,12,11]])
            >>> print(C.front([2]))
            [0 1 3]
            """
            nodes = np.unique(np.asarray(self[startat]))
            front = self.connectedTo(nodes)
            if not add:
                front = np.setdiff1d(front, startat)
            return front

    ######### Creating intermediate levels ###################

        def selectNodes(self, selector):
            """Return a :class:`Connectivity` containing subsets of the nodes.

            Parameters
            ----------
            selector: int array-like
                An object that can be converted to a 1-dim or 2-dim
                int array. Examples are a tuple of local node numbers, or a list
                of such tuples all having the same length.
                Each row of `selector` holds a list of the local node numbers that
                should be retained in the new Connectivity table. As an example,
                if the Connectivity is plex-3 representing triangles, a selector
                [[0,1],[1,2],[2,0]] would extract the edges of the triangle.

            Returns
            -------
            :class:`Connectivity`
                A new Connectivity object with shape
                ``(self.nelems*selector.nelems,selector.nplex)``.
                Duplicate elements created by the selector are retained.
                If the selector has an eltype (for example if it is a
                Connectivity itself), the returned Connectivity will have
                the same eltype.

            Examples
            --------
            >>> Connectivity([[0,1,2],[0,2,1],[0,3,2]]).selectNodes([[0,1],[0,2]])
            Connectivity([[0, 1],
                          [0, 2],
                          [0, 2],
                          [0, 1],
                          [0, 3],
                          [0, 2]])
            """
            sel = Connectivity(selector)
            try:
                eltype = sel.eltype
            except:
                eltype = None
            if sel.size > 0:
                return Connectivity(self[:, sel].reshape(-1, sel.nplex()), eltype=eltype)
            else:
                return Connectivity()

        # BV: should we add a 'unique=False' option to create tables of
        # all intermediate entities without uniquifying?
        #
        def insertLevel(self,selector,permutations='all'):
            """Insert an extra hierarchical level in a Connectivity table.

            A Connectivity table identifies higher hierarchical entities in
            function of lower ones. This method inserts an extra level in the
            hierarchy.
            For example, if you have volumes defined in function of points,
            you can insert an intermediate level of edges, or faces.
            Multiple intermediate level entities may be created from each
            element.

            Parameters
            ----------
            selector: int array-like
                An object that can be converted to a 1-dim or 2-dim
                int array. Examples are a tuple of local node numbers, or a list
                of such tuples all having the same length.
                Each row of `selector` holds a list of the local node numbers that
                should be retained in the new Connectivity table.

                If the Connectivity has an element type, selector can also be a
                single integer specifying one of the hierarchical levels of element
                entities (See the Element class). In that case the selector is
                constructed automatically from self.eltype.getEntities(selector).
            permutations: str
                Defines which permutations of the row data are allowed while still
                considering the rows equal. Equal rows in the intermediate level
                are collapsed into single items. Possible values are:

                - 'none': no permutations are allowed: rows must match the same date
                  at the same positions.
                - 'roll': rolling is allowed. Rows that can be transformed into
                  each other by rolling are considered equal;
                - 'all': any permutation of the same data will be considered an
                  equal row. This is the default.

            Returns
            -------
            hi: :class:`Connectivity`
                A Connecivity defining the original elements
                in function of the intermediate level ones.
            lo: :class:`Connectivity`
                A Connectivity defining the intermediate level
                items in function of the lowest level ones (the original nodes).
                If the ``selector`` has an ``eltype`` attribute, then ``lo`` will
                inherit the same ``eltype`` value.

            The resulting node numbering of the created intermediate entities
            (the `lo` return value) respects the numbering order of the original
            elements and the applied the selector, but in case of collapsing
            duplicate rows, it is undefined which of the collapsed sequences is
            returned.
            Because the precise order of the data in the collapsed rows is lost,
            it is in general not possible to restore the exact original table
            from the two result tables.
            See however :meth:`mesh.Mesh.getBorder` for an application where an
            inverse operation is possible, because the border only contains
            unique rows.

            See also :meth:`mesh.Mesh.combine`, which is an almost inverse operation
            for the general case, if the selector is complete.
            The resulting rows may however be permutations of the original.

            Examples
            --------
            >>> Connectivity([[0,1,2],[0,2,1],[0,3,2]]).\
            insertLevel([[0,1],[1,2],[2,0]])
            (Connectivity([[0, 1, 2],
                          [2, 1, 0],
                          [3, 4, 2]]), Connectivity([[0, 1],
                          [1, 2],
                          [2, 0],
                          [0, 3],
                          [3, 2]]))
            >>> Connectivity([[0,1,2,3]]).\
            insertLevel([[0,1,2],[1,2,3],[0,1,1],[0,0,1],[1,0,0]])
            (Connectivity([[0, 1, 2, 2, 2]]), Connectivity([[0, 1, 2],
                          [1, 2, 3],
                          [0, 1, 1]]))
            """
            from pyformex.elements import elementType

            if at.isInt(selector):
                if hasattr(self, 'eltype'):
                    sel = self.eltype.getEntities(selector)
                else:
                    raise ValueError("Specified an int as selector, but no eltype was defined")
            else:
                sel = Connectivity(selector)
            lo = self.selectNodes(sel)
            if lo.size > 0:
                if sel.testDegenerate().any():
                    LO = lo.copy()
                    # change the double entries to -1
                    LO[np.where(LO[:, :-1] == LO[:, 1:])] = -1
                else:
                    LO = lo
                uniq, uniqid = at.uniqueRowsIndex(LO, permutations=permutations)
                hi = Connectivity(uniqid.reshape(-1, sel.nelems()))
                lo = lo[uniq]
            else:
                hi = lo = Connectivity()
            if hasattr(sel, 'eltype') and sel.eltype is not None:
                lo.eltype = elementType(sel.eltype)
            return hi, lo

        # TODO: This is currently far from general!!!
        # should probably be moved to Mesh/TriSurface if needed there
        def combine(self, lo):
            """Combine two hierarchical Connectivity levels to a single one.

            self and lo are two hierarchical Connectivity tables, representing
            higher and lower level respectively. This means that the elements
            of self hold numbers which point into lo to obtain the lowest level
            items.

            *In the current implementation, the plexitude of lo should be 2!*

            As an example, in a structure of triangles, hi could represent
            triangles defined by 3 edges and lo could represent edges defined
            by 2 vertices. This method will then result in a table
            with plexitude 3 defining the triangles in function of the vertices.

            This is the inverse operation of :meth:`insertLevel` with a selector
            which is complete.
            The algorithm only works if all vertex numbers of an element are
            unique.

            Examples
            --------
            >>> hi,lo = Connectivity([[0,1,2],[0,2,1],[0,3,2]]).\
            insertLevel([[0,1],[1,2],[2,0]])
            >>> hi.combine(lo)
            Connectivity([[0, 1, 2],
                          [0, 2, 1],
                          [0, 3, 2]])

            """
            lo = Connectivity(lo)
            if self.shape[1] < 2 or lo.shape[1] != 2:
                raise ValueError("Can only combine plex>=2 with plex==2")
            elems = lo[self]
            elems1 = np.roll(elems, -1, axis=1)
            for i in range(elems.shape[1]):
                flags = (elems[:, i, 1] != elems1[:, i, 0]) * (elems[:, i, 1] != elems1[:, i, 1])
                elems[flags, i] = np.roll(elems[flags, i], 1, axis=1)
            return Connectivity(elems[:,:, 0])

        def resolve(self):
            """Resolve the connectivity into plex-2 connections.

            Creates a Connectivity table with a plex-2 (edge) connection
            between any two nodes that are connected to a common element.

            There is no point in resolving a plexitude 2 structure.
            Plexitudes lower than 2 can not be resolved.

            Returns a plex-2 Connectivity with all connections between node
            pairs. In each element the nodes are sorted.

            Examples
            --------
            >>> print([ i for i in combinations(range(3),2) ])
            [(0, 1), (0, 2), (1, 2)]
            >>> Connectivity([[0,1,2],[0,2,1],[0,3,2]]).resolve()
            Connectivity([[0, 1],
                          [0, 2],
                          [0, 3],
                          [1, 2],
                          [2, 3]])

            """
            ind = [ i for i in combinations(range(self.nplex()), 2) ]
            hi, lo = self.insertLevel(ind)
            lo.sort(axis=1)
            ind = at.sortByColumns(lo)
            return lo[ind]

        # BV: UNTESTED !
        def reorderNodes(schemes,reverse=False):
            """_Convert Connectivity to/from foreign node numbering schemes.

            The order in which the element's nodes are numbered internally in
            pyFormex may be different than the numbering scheme used in external
            software packages. To allow correct export/import to/from other
            software, the nodes have to be renumbered.
            This function provides such a facility.

            Parameters
            ----------
            schemes: dict
                A dict having pyFormex element names as keys and
                the matching nodal permutation arrays as values. The length of
                the array should match the plexitude of the Connectivity.
            reverse: bool
                If True, the conversion is from external to internal.
                In this case, the Connectivity's eltype is interpreted as the
                pyFormex target element type (and should be set beforehand).

            Returns
            -------
            Connectivity
                If the Connectivity has an element type and `scheme` has a key
                matching the element's name, a Connectivity with the renumbered
                elements is returned.

                - If `reverse` is False (default), the renumbering is done according
                  to the permutation given by the `scheme` value matching the
                  element name and the returned Connectivity will have no element
                  type.

                - If `reverse` is True, the permutation scheme is reversed prior
                  to using it. The target element type is retained in the returned
                  Connectivity.

                If the Connectivity has no element type or `scheme` has no matching
                key, the input Connectivity is returned unchanged.

            """
            if self.eltype is not None:
                eltype = ElementType(self.eltype)
                key = eltype.name()
                if scheme.haskey(key):
                    print('key = %s' % key)
                    trl = scheme[key]
                    print('trl = %s' % trl)
                    elems = self[trl]
                    if not reverse:
                        delattr(self, 'eltype')
                    return elems

            return self

        def sharedNodes(self, elist=None):
            """Return the list of nodes shared by all elements in elist

            Parameters
            ----------
            elist: int array-like
                List of element numbers. If not specified,
                all elements are considered.

            Returns
            -------
            int array
                A 1-dim int array with the list of nodes that are
                common to all elements in the specified list. This array may be
                empty.

            Examples
            --------
            >>> a = Connectivity([[0,1,2],[0,2,1],[0,3,2]])
            >>> a.sharedNodes()
            array([0, 2])
            >>> a.sharedNodes([0,1])
            array([0, 1, 2])

            """
            if elist is None:
                elems = self
            else:
                elems = self[elist]
            m, u = at.multiplicity(elems.ravel())
            return u[m==len(elems)]

        def replic(self, n, inc):
            """Repeat a Connectivity with increasing node numbers.

            Parameters
            ----------
            n: int
                Number of copies to make.
            inc: int
                Increment in node numbers for each copy.

            Returns
            -------
            Connectivity
                A Connectivity with the concatenation of ``n`` replicas of
                ``self``, where the first replica is identical to self and each
                next one has its node numbers increased by ``inc``.

            Examples
            --------
            >>> Connectivity([[0,1,2],[0,2,3]]).replic(2,4)
            Connectivity([[0, 1, 2],
                          [0, 2, 3],
                          [4, 5, 6],
                          [4, 6, 7]])

            """
            return Connectivity(np.concatenate([self+i*inc for i in range(n)]), eltype=self.eltype)

        # TODO: Make this an ElementType method?

        def extrude(self, nnod, degree):
            """Extrude a Connectivity to a higher level Connectivity.

            This method is only available if the Connectivity has an
            element type and the element type is extrudable. See
            :mod:`elements` for more details.

            Parameters
            ----------
            nnod: int
                Node increment for each new node plane. It should be
                higher than the highest node number in self.
            degree: int
                Number of node planes to add. This is also the degree of
                the extrusion. Currently it is limited to 1 or 2.

            The extrusion adds ``degree`` planes of nodes, each with a node
            increment ``nnod``, to the original Connectivity and then selects
            the target nodes from it as defined by the
            ``self.eltype.extruded[degree]`` value.

            Returns
            -------
            Connectivity
                A Connectivity for the extruded element.

            Examples
            --------
            >>> from pyformex.elements import Line2,Line3
            >>> a = Connectivity([[0,1],[1,2]],eltype=Line2).extrude(3,1)
            >>> print(a)
            [[0 1 4 3]
             [1 2 5 4]]
            >>> print(a.eltype.name())
            quad4
            >>> a = Connectivity([[0,1,2],[0,2,3]],eltype=Line3).extrude(4,2)
            >>> print(a)
            [[ 0  2 10  8  1  6  9  4  5]
             [ 0  3 11  8  2  7 10  4  6]]
            >>> print(a.eltype.name())
            quad9

            """
            try:
                eltype, reorder = self.eltype.extruded[degree]
            except:
                try:
                    eltype = self.eltype.name()
                except:
                    eltype = None
                raise ValueError("I don't know how to extrude a Connectivity of eltype '%s' in degree %s" % (eltype, degree))
            # create hypermesh Connectivity
            elems = np.concatenate([self+i*nnod for i in range(degree+1)], axis=-1)
            # Reorder nodes if necessary
            if len(reorder) > 0:
                elems = elems[:, reorder]
            return Connectivity(elems, eltype=eltype)

    #######################################################################
        # class and static methods #

        @staticmethod
        def connect(clist,nodid=None,bias=None,loop=False):
            """Connect nodes from multiple Connectivity objects.

            Parameters
            ----------
            clist: list of Connectivity objects
                The Connectivities to connect.
            nodid: int array-like, optional
                List of node indices, same length as ``clist``. This specifies
                which node of the elements will be used in the connect operation.
            bias: int array-like, optional
                List of element bias values, same length as ``clist``. If provided,
                then element looping will start at this number instead of at zero.
            loop: bool
                If False (default), new element generation will stop as soon as
                the shortest Connectivity runs out of elements. If set to True,
                the shorter lists will wrap around until all elements of all
                Connectivities have been used.

            Returns
            -------
            Connectivity
                A Connectivity with plexitude equal to the number of
                Connectivities in ``clist``. Each element of the new
                Connectivity consist of a node from the corresponding
                element of each of the Connectivities in ``clist``.
                By default this will be the first node of that element,
                but a ``nodid`` list may be given to specify the node index
                to be used for each of the Connectivities.
                Finally, a list of bias values may be given to specify an offset in
                element number for the subsequent Connectivities.
                If loop==False, the length of the Connectivity will be the minimum
                length of the Connectivities in ``clist``, each minus its respective
                bias. If loop=True, the length will be the maximum
                length in of the Connectivities in ``clist``.

            Examples
            --------
            >>> a = Connectivity([[0,1],[2,3],[4,5]])
            >>> b = Connectivity([[10,11,12],[13,14,15]])
            >>> c = Connectivity([[20,21],[22,23]])
            >>> print(Connectivity.connect([a,b,c]))
            [[ 0 10 20]
             [ 2 13 22]]
            >>> print(Connectivity.connect([a,b,c],nodid=[1,0,1]))
            [[ 1 10 21]
             [ 3 13 23]]
            >>> print(Connectivity.connect([a,b,c],bias=[1,0,1]))
            [[ 2 10 22]]
            >>> print(Connectivity.connect([a,b,c],bias=[1,0,1],loop=True))
            [[ 2 10 22]
             [ 4 13 20]
             [ 0 10 22]]
            """
            try:
                m = len(clist)
                for i in range(m):
                    if isinstance(clist[i], Connectivity):
                        pass
                    elif isinstance(clist[i], np.ndarray):
                        clist[i] = Connectivity(clist[i])
                    else:
                        raise TypeError
            except TypeError:
                raise TypeError('Connectivity.connect(): first argument should be a list of Connectivities')

            if not nodid:
                nodid = [ 0 for i in range(m) ]
            if not bias:
                bias = [ 0 for i in range(m) ]
            if loop:
                n = max([ clist[i].nelems() for i in range(m) ])
            else:
                n = min([ clist[i].nelems() - bias[i] for i in range(m) ])
            f = np.zeros((n, m), dtype=at.Int)
            for i, j, k in zip(range(m), nodid, bias):
                v = clist[i][k:k+n, j]
                if loop and k > 0:
                    v = np.concatenate([v, clist[i][:k, j]])
                f[:, i] = np.resize(v, (n))
            return Connectivity(f)

    ############################################################################

    # BV: This could become one of the schemes of Connectivity.reorder
    def _findConnectedLineElems(elems,disconnect=None):
        """Find a single path of connected line elems.

        This function is intended as a helper function for
        :func:`connectedLineElems`. It should probably not be used directly,
        because, as a side-effect, it changes the data in the `elems` argument.
        :func:`connectedLineElems` does not have this inconvenience.

        The function searches a Connectivity table for a chain of elements
        in which the first node of all but the first element is equal to the
        last node of the previous element. To create such a chain, elements may
        be reordered and the node sequence of an element may be reversed.

        Parameters:

        elems: Connectivity-like. Any plexitude is allowed, but only the
          first and the last column are relevant.

        disconnect: optional list/array of integers: node numbers beyond
          which resulting chain can not continue.

        Returns:

        con: a Connectivity with the same shape as the input Connectivity
          `elems`, holding a single chain extracted from the input and filled
          with -1 for the remainder (if any). The chain will not necessarily
          be the longest path. It will however at least contain the first element
          of the input table.

        inv: an int array with two columns and number of rows equal to that of
          `con`. The first column holds the row number in `elems` of the entries
          in `con`. The second column holds a value +1 or -1, flagging whether
          the element is traversed in original direction (+1) in the chain or in
          the reverse direction (-1).

        .. warning:

           As a side-effect, all elements contained in the output Connectivity
           will have their entries in the input table `elems` changed to -1.

        Examples

          >>> con,inv = _findConnectedLineElems([[0,1],[1,2],[0,4],[4,2]])
          >>> print(con)
          [[0 1]
           [1 2]
           [2 4]
           [4 0]]
          >>> print(inv)
          [[ 0  1]
           [ 1  1]
           [ 3 -1]
           [ 2 -1]]

          >>> con,inv = _findConnectedLineElems([[0,1],[1,2],[0,4]])
          >>> print(con)
          [[2 1]
           [1 0]
           [0 4]]
          >>> print(inv)
          [[ 1 -1]
           [ 0 -1]
           [ 2  1]]

          >>> C = Connectivity([[0,1],[0,2],[0,3],[4,5]])
          >>> con,inv = _findConnectedLineElems(C)
          >>> print(con)
          [[ 1  0]
           [ 0  2]
           [-1 -1]
           [-1 -1]]
          >>> print(inv)
          [[ 0 -1]
           [ 1  1]
           [-1  0]
           [-1  0]]
          >>> print(C)
          [[-1 -1]
           [-1 -1]
           [ 0  3]
           [ 4  5]]

        """

        #
        # BV: The side effect in this function could be removed,
        #     because we now have the information on the used elements
        #     in the ind array
        #     That would avoid the need to make a copy in connectedLineElems
        #

        if not isinstance(elems, Connectivity):
            elems = Connectivity(elems)
        #
        # srt will store the sorted connectivity
        #    initialize to -1 to make all invalid
        # ind will store the implicated element numbers and the direction
        #    first column element number, second 1 or -1 (reverse)
        #    first column initialized to -1, second to 0
        # ind is only needed for the return_indices argument in connectedLineElems
        #    but we compute it always for simplicity of the code
        #
        srt = np.zeros_like(elems) - 1
        ind = np.zeros((elems.shape[0], 2), dtype=at.Int)
        ind[:, 0] = -1
        ie = 0
        je = 0
        rev = False
        k = elems[0][0] # remember startpoint
        while True:
            # Store an element that has been found ok
            if rev:
                srt[ie] = elems[je, ::-1]
                ind[ie] = ( je, -1 )
            else:
                srt[ie] = elems[je]
                ind[ie] = ( je, +1 )
            elems[je] = -1 # Done with this one
            j = srt[ie][-1] # remember endpoint
            if j == k or (disconnect is not None and j in disconnect):
                break
            ie += 1

            # Look for the next connected element (only thru first or last node!)
            w = np.where(elems[:, [0, -1]] == j)
            if w[0].size == 0:
                # We've reached the end of a branch
                if disconnect is not None and k in disconnect:
                    # not allowed to revert and continue past start point
                    break
                # Try reversing
                w = np.where(elems[:, [0, -1]] == k)
                if w[0].size == 0:
                    break
                else:
                    j, k = k, j
                    # reverse the table (colums and rows)
                    srt[:ie] = srt[ie-1::-1, ::-1].copy()  # copy needed!!
                    ind[:ie] = ind[ie-1::-1].copy() # rows only
                    ind[:ie, 1] *= -1 # change sign of 2nd column
            je = w[0][0]
            rev = w[-1][0] > 0 #check if the target node is the first or last

        return srt, ind

    # BV: this could become a Connectivity function splitByConnection

    def connectedLineElems(elems,disconnect=None,return_indices=False):
        """Partition a collection of line segments into connected polylines.

        Parameters
        ----------
        elems: intr array with shape (nelems,nplex)
            An 2-dim int array where each row holds the vertex numbers of a
            single line segment. This can e.g.be the ``elems`` Connectivity
            table of a :class:`~mesh.Mesh` with eltype 'line2' (nplex=2).
            But it will work as well for higher order elements like 'line3'
        disconnect: int array-like, optional
            List of vertex numbers where the resulting polylines should be split.
            None of the resulting polylines will have any of the listed vertex
            numbers as an interior vertex. A polyline may start and end
            at such a vertex.
        return_indices: bool
            If True, also return indices identifying the original elements
            contained in each of the resulting polylines.

        Returns
        -------
        conn: list of Connectivity objects
            A list a Connectivity tables of plexitude ``nplex``. The line
            elements in each Connectivity are ordered to form a continuous
            connected segment, i.e. the last vertex of each line element in
            the table is equal to the first vertex of the following element
            (if any).
            The list of connectivity tables is sorted in order of decreasing length.
        ind: list of int arrays, optional
            Only returned if ``return_indices`` is True: a list of tables
            with the same column length as those in ``conn``, and having
            always two columns. The first column contains the original element number
            of the entries, and the second column a value +1 or -1 depending on
            whether the element traversal in the connected segment is in the
            original direction (+1) or the reverse (-1).

        Notes
        -----
        A practical application of the use of the ``disconnect`` paraemter
        is to enforce splitting of the resulting polylines at branching points.
        Branching points can be detected by checking the number of parents.
        ::

            # M is a Mesh of some 'line.' type
            split = where(M.elems.nParents() > 2)[0]
            lines = connectedLineElems(M.elems,disconnect=split)

        Using ``disconnect=None`` would allow polylines to continue over
        the branch points.

        Obviously, from the input elems table and the second return value,
        the first return value could be reconstructed::

          first = [
              where(i[:,-1:] > 0, elems[i[:,0]], elems[i[:,0],::-1])
              for i in second
          ]

        But since the construction of the first list is required by the algorithm,
        it is returned anyway.

        Examples
        --------
        >>> connectedLineElems([[0,1],[1,2],[0,4],[4,2]])
        [Connectivity([[0, 1],
                      [1, 2],
                      [2, 4],
                      [4, 0]])]

        >>> connectedLineElems([[0,1],[1,2],[0,4]])
        [Connectivity([[2, 1],
                      [1, 0],
                      [0, 4]])]

        >>> connectedLineElems([[0,1],[0,2],[0,3],[4,5]])
        [Connectivity([[1, 0],
                      [0, 2]]), Connectivity([[4, 5]]), Connectivity([[0, 3]])]

        >>> connectedLineElems([[0,1],[0,2],[0,3],[4,5]])
        [Connectivity([[1, 0],
                      [0, 2]]), Connectivity([[4, 5]]), Connectivity([[0, 3]])]

        >>> connectedLineElems([[0,1],[0,2],[0,3],[4,5]],return_indices=True)
        ([Connectivity([[1, 0],
                      [0, 2]]), \
        Connectivity([[4, 5]]), Connectivity([[0, 3]])], \
        [array([[ 0, -1],
               [ 1,  1]], dtype=int32), array([[3, 1]], dtype=int32), \
        array([[2, 1]], dtype=int32)])
        >>> connectedLineElems([[0,1,2],[2,0,3],[0,3,1],[4,5,2]])
        [Connectivity([[3, 0, 2],
                      [2, 1, 0],
                      [0, 3, 1]]), Connectivity([[4, 5, 2]])]
        >>> connectedLineElems([[0,1,2],[2,0,3],[0,3,1],[4,5,2]],
        ...     disconnect=[0])
        [Connectivity([[0, 1, 2],
                      [2, 0, 3]]), \
        Connectivity([[4, 5, 2]]), Connectivity([[0, 3, 1]])]
        """
        elems = Connectivity(elems).copy() # make copy to avoid side effects
        elnrs =  np.arange(elems.shape[0]) # needed to return indices
        parts = []
        chains = []
        while elems.size != 0:
            loop, ind = _findConnectedLineElems(elems,disconnect)
            ind[:, 0] = elnrs[ind[:, 0]]
            parts.append(loop[(loop!=-1).any(axis=1)])
            chains.append(ind[(loop!=-1).any(axis=1)])
            todo = (elems!=-1).any(axis=1)
            elems = elems[todo]
            elnrs = elnrs[todo]
        # sort according to decreasing number of elements
        nel = [ p.nelems() for p in parts ]
        srt = np.argsort(nel)[::-1]
        parts = [ parts[i] for i in srt ]
        chains = [ chains[i] for i in srt ]
        if return_indices:
            return parts, chains
        else:
            return parts

    # Private functions for adjacency multiprocessing

    def _elem_adj(inv,els,check):
        """Return elem adj for (part of) the elements"""
        adj = inv[els].reshape((els.shape[0], -1))
        return Adjacency(adj,check_max=check)

    def _elem_adj_multi(inv,els,nproc=-1):
        from pyformex import multi
        datablocks = at.splitar(els, nproc)
        datalen = [0] + [d.shape[0] for d in datablocks]
        shift = np.array(datalen[:-1]).cumsum()
        tasks = [(_elem_adj, (inv, e, False)) for e,s in zip (datablocks,shift)]
        return multi.multitask(tasks, nproc=nproc)

    # End