def nodes_not_in_edge(nodes, edge):
    """Given a list of nodes numbers, return the indices of those which do
    not belong to the specified edge

    Parameters
    ----------
    nodes : array or list
        The list of node numbers
    edge: array or list
        The edge, defined by two (or more) node number

    Return
    ------
    non_member_nodes : list
        The nodes which are not a member of the given edge
    """
    return [node_index for node_index, node_num in enumerate(nodes)
            if node_num not in edge] def shared_nodes(nodes1, nodes2):
    """Return all the nodes shared by two polyhedra

    Parameters
    ----------
    nodes1 : array/list
        A list of node numbers
    nodes2 : array/list
        A list of node numbers

    Returns
    -------
    common_nodes : list
        A list of nodes found in both lists
    """
    return [node for node in nodes1 if node in nodes2] class cached_property(object):
    """
    A property that is only computed once per instance and then replaces itself
    with an ordinary attribute. Deleting the attribute resets the property.

    Taken from https://github.com/pydanny/cached-property/blob/master/cached_property.py
    Original source: https://github.com/bottlepy/bottle/commit/fa7733e075da0d790d809aa3d2f53071897e6f76
    Copyright under MIT License
    """

    def __init__(self, func):
        self.__doc__ = getattr(func, '__doc__')
        self.func = func

    def __get__(self, obj, cls):
        if obj is None:
            return self
        value = obj.__dict__[self.func.__name__] = self.func(obj)
        return value def inc_slice(s, inc):
    """Increment a slice so that it starts at the current stop, and the current
    stop is incremented by some amount"""
    return slice(s.stop, s.stop+inc) class Identified(object):
    """An object which can be uniquely identified by an id number. It is
    assumed that any object which subclasses Identified is immutable, so that
    its id can be used for caching complex results which depend on this object.
    """

    def __init__(self):
        self.id = uuid.uuid4()

    def __hash__(self):
        return self.id.__hash__()

    def __eq__(self, other):
        return hasattr(other, 'id') and (self.id == other.id)

    def __repr__(self):
        "Represent the object by its id, in addition to its memory address"
        return ("<%s at 0x%08x with id %s>" % (str(self.__class__)[8:-2],
                                               id(self),
                                               str(self.id))) def memoize(obj):
    """A simple decorator to memoize function calls. Pays particular attention
    to numpy arrays and objects which are subclasses of Identified. It is
    assumed that in such cases, the object does not change if its `id` is the
    same"""
    cache = obj.cache = {}

    def get_key(item):
        if isinstance(item, (six.string_types, numbers.Number)):
            return item
        elif isinstance(item, Identified):
            return str(item.id)
        elif isinstance(item, np.ndarray):
            return item.tostring()
        else:
            return str(item)

    @functools.wraps(obj)
    def memoizer(*args, **kwargs):
        key_arg = tuple(get_key(arg) for arg in args)
        key_kwarg = tuple((kw, get_key(arg)) for (kw, arg)
                          in kwargs.items())
        key = (key_arg, key_kwarg)

        if key not in cache:
            cache[key] = obj(*args, **kwargs)
        return cache[key]
    return memoizer def equivalence(relations):
    """Determine the equivalence classes between objects

    Following numerical recipes section 8.6

    Parameters
    ----------
    relations: list
        Each element of the list is a tuple containing the identities of two
        equivalent items. Each item can be any hashable type

    Returns
    -------
    class_items: list of set
        Each set
    """

    # first put each item in its own equivalence class
    classes = {}
    for j, k in relations:
        classes[j] = j
        classes[k] = k

    for relation in relations:
        j, k = relation

        # track the anscestor of each
        while classes[j] != j:
            j = classes[j]

        while classes[k] != k:
            k = classes[k]

        # if not already related, then relate items
        if j != k:
            classes[j] = k

    # The final sweep
    for j in classes.keys():
        while classes[j] != classes[classes[j]]:
            classes[j] = classes[classes[j]]

    # Now reverse the storage arrangement, so that all items of the same
    # class are grouped together into a set
    classes_reverse = defaultdict(set)
    for item, item_class in classes.items():
        classes_reverse[item_class].add(item)

    # the class names are arbitrary, so just return the list of sets
    return list(classes_reverse.values()) class MeshError(Exception):
    "An exeception indicating a failure generating or reading the mesh"
    pass class IntegrationRule(Identified):
class DunavantRule(IntegrationRule):
class GaussLegendreRule(IntegrationRule):
class TrapezoidalRule(IntegrationRule):
class Contour(object):
class CircularContour(Contour):
class RectangularContour(Contour):
class ExternalModeContour(Contour):
class EllipticalContour(Contour):
    def __len__(self):
    def __repr__(self):
    def __iter__(self):
    def __init__(self, order):
    def __init__(self, order):
    def __init__(self, order):
def cartesian_to_barycentric(r, nodes):
def triangle_electric_dipole(vertices, xi_eta, weights):
def sphere_fibonacci(num_points, cartesian=False):
    def points_inside(self, points):
    def __init__(self, centre, radius, integration_rule=TrapezoidalRule(20)):
    def __iter__(self):
    def __len__(self):
    def __init__(self, s_min, s_max, integration_rule=GaussLegendreRule(20)):
    def __iter__(self):
    def __len__(self):
    def __init__(self, corner, integration_rule=GaussLegendreRule(20),
                 overlap_axes=None, avoid_origin=None):
    def __iter__(self):
    def __len__(self):
    def __init__(self, radius_real, radius_imag, offset_real, offset_imag,
                 integration_rule=GaussLegendreRule(20)):
    def __iter__(self):
    def __len__(self):
    T = np.array(((nodes[0, 0] - nodes[2, 0], nodes[1, 0] - nodes[2, 0]),
                 (nodes[0, 1] - nodes[2, 1], nodes[1, 1] - nodes[2, 1]))) class Part(Identified):
    """A part which has been placed into the simulation, and which can be
    modified"""

    def __init__(self, location=None):
        super(Part, self).__init__()

        self.initial_location = location
        self.transformation_matrix = np.empty((4, 4))
        self.reset()
        self.parent_ref = None

    def reset(self):
        """Reset this part to the default values of the original `Mesh`
        from which this `Part` was created, and move it to its initial
        location (the origin if not specified).
        """
        self.transformation_matrix[:] = np.eye(4)
        if self.initial_location is not None:
            self.translate(self.initial_location)
        self._position_updated()

    def _position_updated(self):
        "Called when position is updated"
        pass

    @property
    def complete_transformation(self):
        """The complete transformation matrix, which takes into account
        the transformation matrix of all parents"""
        if self.parent_ref is None:
            return self.transformation_matrix
        else:
            return self.parent_ref().complete_transformation.dot(
                                                    self.transformation_matrix)

    def translate(self, offset_vector):
        """Translate a part by an arbitrary offset vector

        Care needs to be take if this puts an object in a different layer
        """
        translation = np.eye(4)
        translation[:3, 3] = offset_vector

        self.transformation_matrix[:] = translation.dot(
                                                    self.transformation_matrix)
        self._position_updated()

    def rotate(self, axis, angle):
        """
        Rotate about an arbitrary axis

        Parameters
        ----------
        axis : ndarray
            the vector about which to rotate
        angle : number
            angle of rotation in degrees

        Algorithm taken from
        http://en.wikipedia.org/wiki/Euler%E2%80%93Rodrigues_parameters
        """

        # TODO: enable rotation about arbitrary coordinates, and about the
        # centre of the object

        axis = np.array(axis, np.float64)
        axis /= np.sqrt(np.dot(axis, axis))

        angle = angle*np.pi/180.0

        a = np.cos(0.5*angle)
        b, c, d = axis*np.sin(0.5*angle)

        matrix = np.array([[a**2 + b**2 - c**2 - d**2,
                            2*(b*c - a*d), 2*(b*d + a*c), 0],
                           [2*(b*c + a*d), a**2 + c**2 - b**2 - d**2,
                            2*(c*d - a*b), 0],
                           [2*(b*d - a*c), 2*(c*d + a*b),
                            a**2 + d**2 - b**2 - c**2, 0],
                           [0, 0, 0, 1.0]])

        self.transformation_matrix[:] = matrix.dot(self.transformation_matrix)
        self._position_updated()
import logging
import numpy as np
from collections import namedtuple
from scipy.sparse import lil_matrix
from openmodes.mesh import nodes_not_in_edge, shared_nodes
from openmodes.helpers import (cached_property, inc_slice, Identified, memoize,
                               equivalence, MeshError)
from openmodes.integration import triangle_centres
from openmodes.external.ordered_set import OrderedSet
from openmodes.parts import Part
        tri_m = np.empty(num_basis, np.int32)
        # internal index of free node of T+
        node_p = np.empty(num_basis, np.int32)
        # internal index of free node of T-
        node_m = np.empty(num_basis, np.int32)

        for basis_count, edge_count in enumerate(shared_edge_indices):
            # set the RWG basis function triangles
            tri_p[basis_count] = triangles_shared_by_edges[edge_count][0]
            tri_m[basis_count] = triangles_shared_by_edges[edge_count][1]

            # determine the indices of the unshared nodes, indexed within the
            # sharing triangles (i.e. 0, 1 or 2)
            node_p[basis_count] = nodes_not_in_edge(
                       mesh.polygons[tri_p[basis_count]], edges[edge_count])[0]
            node_m[basis_count] = nodes_not_in_edge(
                       mesh.polygons[tri_m[basis_count]], edges[edge_count])[0]

        self.rwg = RWG(tri_p, tri_m, node_p, node_m)
        self.sections = (num_basis,)

        logging.info("Constructing %d RWG basis functions over %d faces"
                     % (num_basis, len(mesh.polygons)))

    @cached_property
    def transformation_matrices(self):
        """Returns the sparse transformation matrix to turn quantities
        defined on faces to loop-star basis

        For vector quantities, assumes that the face-based quantity has been
        packed to a 2D array of size (n_basis*3, n_basis*3)

        For scalar quantities, assumes that the face-based quantity has been
        packed to a 2D array of size (n_basis, n_basis)
        """

        num_basis = len(self)
        num_tri = len(self.mesh.polygons)
        # scalar_transform = np.zeros((num_basis, num_tri), np.float64)
        # vector_transform=np.zeros((num_basis, 3*num_tri), np.float64)
        scalar_transform = lil_matrix((num_basis, num_tri))
        vector_transform = lil_matrix((num_basis, 3*num_tri))

        for basis_count, (tri_p, tri_m, node_p, node_m) in enumerate(self):
            scalar_transform[basis_count, tri_p] = 1.0
            scalar_transform[basis_count, tri_m] = -1.0

            vector_transform[basis_count, tri_p*3+node_p] = 1.0
            vector_transform[basis_count, tri_m*3+node_m] = -1.0

        return vector_transform.tocsr(), scalar_transform.tocsr()

    def __len__(self):
        return len(self.rwg.tri_p)

    def __getitem__(self, index):
        return RWG(self.rwg.tri_p[index], self.rwg.tri_m[index],
                   self.rwg.node_p[index], self.rwg.node_m[index])


def construct_stars(mesh, edges, triangles_shared_by_edges, sharing_count):
    """Construct star basis functions on a triangular mesh. The star
    corresponding to one triangle faces will be arbitrarily dropped"""

    num_tri = len(mesh.polygons)

    shared_edge_indices = np.where(sharing_count == 2)[0]

    tri_p = [list() for _ in range(num_tri)]
    tri_m = [list() for _ in range(num_tri)]
    node_p = [list() for _ in range(num_tri)]
    node_m = [list() for _ in range(num_tri)]

    # Go through shared edges, and update both star-basis functions
    # to add the influence of this shared edge.
    for edge_count in shared_edge_indices:
        tri1, tri2 = triangles_shared_by_edges[edge_count]

        tri_p[tri1].append(tri1)
        tri_p[tri2].append(tri2)

        tri_m[tri1].append(tri2)
        tri_m[tri2].append(tri1)

        node1 = nodes_not_in_edge(mesh.polygons[tri1],
                                  edges[edge_count])[0]
        node2 = nodes_not_in_edge(mesh.polygons[tri2],
                                  edges[edge_count])[0]
        node_p[tri1].append(node1)
        node_p[tri2].append(node2)

        node_m[tri1].append(node2)
        node_m[tri2].append(node1)

    return RWG(tri_p[:-1], tri_m[:-1], node_p[:-1], node_m[:-1])


def construct_loop(loop_triangles, polygons):
    """Construct a single loop basis function corresponding to a single inner
    node of a triangular mesh"""

    tri_p = []
    tri_m = []
    node_p = []
    node_m = []

    # The first triangle to be picked doesn't get removed from
    # the list, as we will need to refer to it again
    current_triangle = loop_triangles.pop()
    first_triangle = current_triangle

    while len(loop_triangles) > 0:
        for triangle_count, next_triangle in enumerate(loop_triangles):
            shared = shared_nodes(polygons[next_triangle],
                                  polygons[current_triangle])
            if len(shared) == 2:
                break
        loop_triangles.pop(triangle_count)

        # find the unshared nodes

        free_current = nodes_not_in_edge(polygons[current_triangle], shared)[0]