# Source code for trisurface

    #
    ##
    ##  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/.
    ##
    """Operations on triangulated surfaces.

    A triangulated surface is a surface consisting solely of triangles.
    Any surface in space, no matter how complex, can be approximated with
    a triangulated surface.
    """
    from __future__ import absolute_import, division, print_function

    import os
    import numpy as np

    import pyformex as pf
    from pyformex import fileread, filewrite, geomtools, inertia, utils
    from pyformex.coords import Coords
    from pyformex.connectivity import Connectivity, connectedLineElems
    from pyformex.mesh import Mesh
    from pyformex.formex import Formex
    from pyformex.geometry import Geometry
    from pyformex.arraytools import *

    #
    # gts commands used:
    #   in Debian package: stl2gts gts2stl gtscheck
    #   not in Debian package: gtssplit gtscoarsen gtsrefine gtssmooth gtsinside
    #

    # Conversion of surface file formats

    #
    # BV: the following functions have to be checked for their need
    # and opportunity, and replaced by more general infrastrucuture
    #
    def adjacencyArrays(elems,nsteps=1):
        """Create adjacency arrays for 2-node elements.

        elems is a (nr,2) shaped integer array.
        The result is a list of adjacency arrays, where row i of adjacency array j
        holds a sorted list of the nodes that are connected to node i via a shortest
        path of j elements, padded with -1 values to create an equal list length
        for all nodes.
        This is: [adj0, adj1, ..., adjj, ... , adjn] with n=nsteps.

        Examples

        >>> for a in adjacencyArrays([[0,1],[1,2],[2,3],[3,4],[4,0]],3):
        ...     print(a)
        [[0]
         [1]
         [2]
         [3]
         [4]]
        [[1 4]
         [0 2]
         [1 3]
         [2 4]
         [0 3]]
        [[2 3]
         [3 4]
         [0 4]
         [0 1]
         [1 2]]
        []

        """
        #utils.warn("depr_adjacencyArrays")
        from pyformex import adjacency
        elems = Connectivity(elems)
        if len(elems.shape) != 2 or elems.shape[1] != 2:
            raise ValueError("""Expected a set of 2-node elements.""")
        if nsteps < 1:
            raise ValueError("""The shortest path should be at least 1.""")
        # Construct table of nodes connected to each node
        adj1 = elems.adjacency('n')
        m = adj1.shape[0]
        adj = [ np.arange(m).reshape(-1, 1), adj1 ]
        nodes = adj1
        step = 2
        while step <= nsteps and nodes.size > 0:
            # Determine adjacent nodes
            t = nodes < 0
            nodes = adj1[nodes]
            nodes[t] = -1
            nodes = nodes.reshape((m, -1))
            nodes = adjacency._reduceAdjacency(nodes)
            # Remove nodes of lower adjacency
            ladj = np.concatenate(adj[-2:], -1)
            t = [ np.in1d(n, l, assume_unique=True) for n, l in zip (nodes, ladj) ]
            t = np.asarray(t)
            nodes[t] = -1
            nodes = adjacency._sortAdjacency(nodes)
            # Store current nodes
            adj.append(nodes)
            step += 1
        return adj

    #
    # TODO: this should unzip compressed input files and zip compressed output
    #
    def stlConvert(stlname,outname=None,binary=False,options='-d'):
        """Transform an .stl file to .off or .gts or binary .stl format.

        Parameters:

        - `stlname`: name of an existing .stl file (either ascii or binary).
        - `outname`: name of the output file. The extension defines the format
          and should be one of '.off', '.gts', '.stl', '.stla', or .stlb'.
          As a convenience, if a file extension only is given (other than '.stl'),
          then the outname will be constructed by changing the extension of the
          input `stlname`.
        - `binary`: if the extension of `outname` is '.stl', defines whether
          the output format is a binary or ascii STL format.

        If the outname file exists and its mtime is more recent than the stlname,
        the outname file is considered uptodate and the conversion program will
        not be run.

        The conversion program will be choosen depending on the extension.
        This uses the external commands 'admesh' or 'stl2gts'.

        The return value is a tuple of the output file name, the conversion
        program exit code (0 if succesful) and the stdout of the conversion
        program (or a 'file is already uptodate' message).
        """
        if not outname:
            outname = pf.cfg.get('surface/stlread', '.off')
        if outname.startswith('.'):
            outname = utils.changeExt(stlname, outname)
        if os.path.exists(outname) and utils.mtime(stlname) < utils.mtime(outname):
            return outname, 0, "File '%s' seems to be up to date" % outname

        ftype = utils.fileTypeFromExt(outname)
        if ftype == 'stl' and binary:
            ftype = 'stlb'

        if ftype == 'off':
            utils.hasExternal('admesh')
            cmd = "admesh %s --write-off '%s' '%s'" % (options, outname, stlname)
        elif ftype in [ 'stl', 'stla' ]:
            utils.hasExternal('admesh')
            cmd = "admesh %s -a '%s' '%s'" % (options, outname, stlname)
        elif ftype == 'stlb':
            utils.hasExternal('admesh')
            cmd = "admesh %s -b '%s' '%s'" % (options, outname, stlname)
        elif ftype == 'gts':
            cmd = "stl2gts < '%s' > '%s'" % (stlname, outname)
        else:
            return outname, 1, "Can not convert file '%s' to '%s'" % (stlname, outname)

        P = utils.command(cmd, shell=True)
        return outname, P.sta, P.out

    # Input of surface file formats

    def read_stl(fn,intermediate=None):
        """Read a surface from .stl file.

        This is done by first coverting the .stl to .gts or .off format.
        The name of the intermediate file may be specified. If not, it will be
        generated by changing the extension of fn to '.gts' or '.off' depending
        on the setting of the 'surface/stlread' config setting.

        Return a coords,edges,faces or a coords,elems tuple, depending on the
        intermediate format.
        """
        ofn, sta, out = stlConvert(fn, intermediate)
        if sta:
            pf.debug("Error during conversion of file '%s' to '%s'" % (fn, ofn))
            pf.debug(out)
            return ()

        ftype,compr = utils.fileTypeComprFromExt(ofn)
        if ftype == 'off':
            return fileread.read_off(ofn)
        elif ftype == 'gts':
            return fileread.read_gts(ofn)

    def curvature(coords,elems,edges,neighbours=1):
        """Calculate curvature parameters at the nodes.

        Algorithms based on Dong and Wang 2005; Koenderink and Van Doorn 1992.
        This uses the nodes that are connected to the node via a shortest
        path of 'neighbours' edges.
        Eight values are returned: the Gaussian and mean curvature, the
        shape index, the curvedness, the principal curvatures and the
        principal directions.
        """
        # calculate n-ring neighbourhood of the nodes (n=neighbours)
        adj = adjacencyArrays(edges, nsteps=neighbours)[-1]
        adjNotOk = adj<0
        # for nodes that have less than three adjacent nodes, remove the adjacencies
        adjNotOk[(adj>=0).sum(-1) <= 2] = True
        # calculate unit length average normals at the nodes p
        # a weight 1/|gi-p| could be used (gi=center of the face fi)
        p = coords
        n = geomtools.averageNormals(coords, elems, atNodes=True)
        # double-precision: this will allow us to check the sign of the angles
        p = p.astype(float64)
        n = n.astype(float64)
        vp = p[adj] - p[:, newaxis]
        vn = n[adj] - n[:, newaxis]
        # where adjNotOk, set vectors = [0.,0.,0.]
        # this will result in NaN values
        vp[adjNotOk] = 0.
        vn[adjNotOk] = 0.
        # calculate unit length projection of vp onto the tangent plane
        t = geomtools.projectionVOP(vp, n[:, newaxis])
        t = normalize(t)
        # calculate normal curvature
        k = dotpr(vp, vn)/dotpr(vp, vp)
        # calculate maximum normal curvature and corresponding coordinate system
        try:
            imax = nanargmax(k, -1)
            kmax =  k[arange(len(k)), imax]
            tmax = t[arange(len(k)), imax]
        except: # bug with numpy.nanargmax: cannot convert float NaN to integer
            kmax = resize(NaN, (k.shape[0]))
            tmax = resize(NaN, (t.shape[0], 3))
            w = ~(isnan(k).all(1))
            imax = nanargmax(k[w], -1)
            kmax[w] =  k[w, imax]
            tmax[w] =  t[w, imax]
        tmax1 = tmax
        tmax2 = cross(n, tmax1)
        tmax2 = normalize(tmax2)
        # calculate angles (tmax1,t)
        theta, rot = geomtools.rotationAngle(repeat(tmax1[:, newaxis], t.shape[1], 1), t, angle_spec=RAD)
        theta = theta.reshape(t.shape[:2])
        rot = rot.reshape(t.shape)
        # check the sign of the angles
        d =  dotpr(rot, n[:, newaxis])/(length(rot)*length(n)[:, newaxis]) # divide by length for round-off errors
        cw = np.isclose(d, [-1.])
        theta[cw] = -theta[cw]
        # calculate coefficients
        a = kmax
        a11 = nansum(cos(theta)**2*sin(theta)**2, -1)
        a12 = nansum(cos(theta)*sin(theta)**3, -1)
        #a21 = a12
        a22 = nansum(sin(theta)**4, -1)
        a13 = nansum((k-a[:, newaxis]*cos(theta)**2)*cos(theta)*sin(theta), -1)
        a23 = nansum((k-a[:, newaxis]*cos(theta)**2)*sin(theta)**2, -1)
        denom = (a11*a22-a12**2)
        b = (a13*a22-a23*a12)/denom
        c = (a11*a23-a12*a13)/denom
        # calculate the Gaussian and mean curvature
        K = a*c-b**2/4
        H = (a+c)/2
        # calculate the principal curvatures and principal directions
        k1 = H+sqrt(H**2-K)
        k2 = H-sqrt(H**2-K)
        theta0 = 0.5*arcsin(b/(k2-k1))
        w = apply_along_axis(np.isclose, 0, -b, 2*(k2-k1)*cos(theta0)*sin(theta0))
        theta0[w] = pi-theta0[w]
        e1 = cos(theta0)[:, newaxis]*tmax1+sin(theta0)[:, newaxis]*tmax2
        e2 = cos(theta0)[:, newaxis]*tmax2-sin(theta0)[:, newaxis]*tmax1
        # calculate the shape index and curvedness
        S = 2./pi*arctan((k1+k2)/(k1-k2))
        C = square((k1**2+k2**2)/2)
        return K, H, S, C, k1, k2, e1, e2

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

    def fillBorder(border,method='radial',dir=None):
        """Create a surface inside a given closed border line.

        The border line is a closed polygonal line and can be specified as
        one of the following:

        - a closed PolyLine,
        - a 2-plex Mesh, with a Connectivity table such that the elements
          in order form a closed polyline,
        - a simple Coords specifying the subsequent vertices of the polygonal
          border line.

        The return value is a TriSurface filling the hole inside the border.

        There are currently two fill methods available:

        - 'radial': this method adds a central point and connects all border
          segments with the center to create triangles.
        - 'border': this method creates subsequent triangles by connecting the
          endpoints of two consecutive border segments and thus works its way
          inwards until the hole is closed. Triangles are created at the line
          segments that form the smallest angle.

        The 'radial' method produces nice results if the border is relative smooth,
        nearly convex and nearly planar. It adds an extra point though, which may
        be unwanted. On irregular 3D borders there is a high change that the
        result contains intersecting triangles.

        This 'border' method is slower on large borders, does not introduce any
        new point and has a better chance of avoiding intersecting triangles
        on irregular 3D borders.

        The resulting surface can be checked for intersecting triangles by the
        :meth:`check` method.

        .. note :: Because the 'border' does not create any new points, the
          returned surface will use the same point coordinate array as the input
          object.
        """
        from pyformex.plugins.curve import PolyLine
        if isinstance(border, Mesh) and border.nplex()==2:
            if method == 'radial':
                border = border.compact()
            coords = border.coords
            elems = border.elems[:, 0]
        elif isinstance(border, PolyLine):
            coords = border.coords
            elems = None
        elif isinstance(border, Coords):
            coords = border.reshape(-1, 3)
            elems = None
        else:
            raise ValueError("Expected a 2-plex Mesh, a PolyLine or a Coords array as first argument")

        if elems is None:
            elems = arange(coords.shape[0])

        n = elems.shape[0]
        if n < 3:
            raise ValueError("Expected at least 3 points.")

        if method == 'radial':
            coords = Coords.concatenate([coords, coords.center()])
            elems = column_stack([elems, roll(elems, -1), n*ones(elems.shape[0], dtype=Int)])

        elif method == 'border':
            # creating elems array at once (more efficient than appending)
            tri = -ones((n-2, 3), dtype=Int)
            # compute all internal angles
            x = coords[elems]
            e = arange(n)
            v = roll(x, -1, axis=0) - x
            v = normalize(v)
            c = vectorPairCosAngle(roll(v, 1, axis=0), v)
            # loop in order of smallest angles
            itri = 0
            while n > 3:
                # find minimal angle
                j = c.argmin()
                i = (j - 1) % n
                k = (j + 1) % n
                tri[itri] = [ e[i], e[j], e[k]]
                # remove the point j of triangle i,j,k
                # recompute adjacent angles of edge i,k
                ii = (i-1) % n
                v1 = normalize([ v[e[ii]], x[e[k]] - x[e[i]] ])
                v2 = normalize([ x[e[k]] - x[e[i]], v[e[k]] ])
                cnew = vectorPairCosAngle(v1, v2)
                c = roll(concatenate([cnew, roll(c, 1-j)[3:]]), j-1)
                e = roll(roll(e, -j)[1:], j)
                n -= 1
                itri += 1
            tri[itri] = e
            elems = elems[tri]

        elif method == 'planar':
            from pyformex.plugins import polygon as pg
            x = coords[elems]
            e = arange(x.shape[0])

            if dir is None:
                dir = geomtools.smallestDirection(x)

            X, C, A, a = pg.projected(x, dir)
            P = pg.Polygon(X)
            if P.area() < 0.0:
                P = P.reverse()
                e = reverseAxis(e, 0)
            S = P.fill()
            e = e[S.elems]
            elems = elems[e]

        else:
            raise ValueError("Strategy should be either 'radial', 'border' or 'planar'")

        return TriSurface(coords, elems)

    ############################################################################
    # The TriSurface class

    class TriSurface(Mesh):
        """A class representing a triangulated 3D surface.

        The surface contains `ntri` triangles, each having 3 vertices with
        3 coordinates. The surface can be initialized from one of the following:

        - a (ntri,3,3) shaped array of floats
        - a Formex with plexitude 3
        - a Mesh with plexitude 3
        - an (ncoords,3) float array of vertex coordinates and
          an (ntri,3) integer array of vertex numbers
        - an (ncoords,3) float array of vertex coordinates,
          an (nedges,2) integer array of vertex numbers,
          an (ntri,3) integer array of edges numbers.

        Additionally, a keyword argument prop= may be specified to
        set property values.
        """

        _exclude_members_ = [ 'intersectionWithLines' ]

        def __init__(self,*args,**kargs):
            """Create a new surface."""
            self._areas = self._fnormals = None
            self.adj = None
            if hasattr(self, 'edglen'):
                del self.edglen

            if len(args) == 0:
                Mesh.__init__(self, [], [], None, 'tri3')
                return  # an empty surface

            if len(args) == 1:
                # argument should be a suitably structured geometry object
                # TriSurface, Mesh, Formex, Coords, ndarray, ...
                a = args[0]

                if isinstance(a, Mesh):
                    if a.nplex() != 3 or a.elName() != 'tri3':
                        raise ValueError("Only meshes with plexitude 3 and eltype 'tri3' can be converted to TriSurface!")
                    Mesh.__init__(self, a.coords, a.elems, a.prop, 'tri3')

                else:
                    if not isinstance(a, Formex):
                        # something that can be converted to a Formex
                        try:
                            a = Formex(a)
                        except:
                            raise ValueError("Can not convert objects of type %s to TriSurface!" % type(a))

                    # now a is a Formex
                    if a.nplex() != 3:
                        raise ValueError("Expected an object with plexitude 3!")

                    coords, elems = a.coords.fuse()
                    Mesh.__init__(self, coords, elems, a.prop, 'tri3')

            else:
                # arguments are (coords,elems) or (coords,edges,faces)
                coords = Coords(args[0])
                if len(coords.shape) != 2:
                    raise ValueError("Expected a 2-dim coordinates array")

                if len(args) == 2:
                    # arguments are (coords,elems)
                    elems = Connectivity(args[1], nplex=3)
                    Mesh.__init__(self, coords, elems, None, 'tri3')

                elif len(args) == 3:
                    # arguments are (coords,edges,faces)
                    edges = Connectivity(args[1], nplex=2)

                    if edges.size > 0 and edges.max() >= coords.shape[0]:
                        raise ValueError("Some vertex number is too high")

                    faces = Connectivity(args[2], nplex=3)

                    if faces.max() >= edges.shape[0]:
                        raise ValueError("Some edge number is too high")

                    elems = faces.combine(edges)
                    Mesh.__init__(self, coords, elems, None, 'tri3')

                    # since we have the extra data available, keep them
                    self.edges = edges
                    self.elem_edges = faces

                else:
                    raise RuntimeError("Too many positional arguments")

            if 'prop' in kargs:
                self.setProp(kargs['prop'])

        def __setstate__(self, state):
            """Set the object from serialized state.

            This allows to read back old pyFormex Project files where the
            Surface class did not set an element type.
            """
            if 'areas' in state:
                state['_areas'] = state['areas']
                del state['areas']
            self.__dict__.update(state)
            self.setType('tri3')

    ###########################################################################
        #
        #   Return information about a TriSurface
        #

        def nedges(self):
            """Return the number of edges of the TriSurface."""
            return self.getEdges().shape[0]

        def nfaces(self):
            """Return the number of faces of the TriSurface."""
            return self.getElemEdges().shape[0]

        def vertices(self):
            """Return the coordinates of the nodes of the TriSurface."""
            return self.coords

        def shape(self):
            """Return the number of points, edges, faces of the TriSurface."""
            return self.ncoords(), self.nedges(), self.nfaces()

        def getElemEdges(self):
            """Get the faces' edge numbers."""
            if self.elem_edges is None:
                self.elem_edges, self.edges = self.elems.insertLevel(1)
            return self.elem_edges

    ###########################################################################
        #
        #   Operations that change the TriSurface itself
        #
        #  Make sure that you know what you're doing if you use these
        #
        #
        # Changes to the geometry should by preference be done through the
        # __init__ function, to ensure consistency of the data.
        # Convenience functions are defined to change some of the data.
        #

        def setCoords(self, coords):
            """Change the coords."""
            self.__init__(coords, self.elems, prop=self.prop)
            return self

        def setElems(self, elems):
            """Change the elems."""
            self.__init__(self.coords, elems, prop=self.prop)

        def setEdgesAndFaces(self, edges, faces):
            """Change the edges and faces."""
            self.__init__(self.coords, edges, faces, prop=self.prop)

        def append(self, S):
            """Merge another surface with self.

            This just merges the data sets, and does not check
            whether the surfaces intersect or are connected!
            This is intended mostly for use inside higher level functions.
            """
            coords = concatenate([self.coords, S.coords])
            elems = concatenate([self.elems, S.elems+self.ncoords()])
            ## What to do if one of the surfaces has properties, the other one not?
            ## The current policy is to use zero property values for the Surface
            ## without props
            prop = None
            if self.prop is not None or S.prop is not None:
                if self.prop is None:
                    self.prop = zeros(shape=self.nelems(), dtype=Int)
                if S.prop is None:
                    p = zeros(shape=S.nelems(), dtype=Int)
                else:
                    p = S.prop
                prop = concatenate((self.prop, p))
            self.__init__(coords, elems, prop=prop)

    ###########################################################################
        #
        #   read and write
        #

        @classmethod
        def read(clas,fn,ftype=None):
            """Read a surface from file.

            If no file type is specified, it is derived from the filename
            extension.
            Currently supported file types:

            - .off
            - .gts
            - .stl (ASCII or BINARY)
            - .neu (Gambit Neutral)
            - .smesh (Tetgen)

            Compressed (gzip or bzip2) files are also supported. Their names
            should be the normal filename with '.gz' or '.bz2' appended.
            These files are uncompressed on the fly during the reading and
            the uncompressed versions are deleted after reading.

            The file type can be specified explicitely to handle file names
            where the extension does not directly specify the file type.
            """
            if ftype is None:
                ftype,compr = utils.fileTypeComprFromExt(fn)
            else:
                ftype,compr = utils.fileTypeComprFromExt('a.'+ftype)

            if ftype == 'off':
                data = fileread.read_off(fn)
            elif ftype == 'gts':
                data = fileread.read_gts(fn)
            elif ftype == 'stl':
                try:
                    # first try binary format
                    data, color = fileread.read_stl_bin(fn)
                    # data has 4 items per triangle: normal + 3 vertices
                    S = TriSurface(data[:, 1:])
                    if color:
                        S.attrib(color = color)
                    return S
                except:
                    print("Could not read as binary stl, will try conversion")
                    data = read_stl(fn)
            # The remainder should probably disappear
            # The following do not support compression yet
            elif ftype == 'smesh' and not compr:
                from pyformex.plugins import tetgen
                data = tetgen.readSurface(fn)
            elif ftype == 'neu' and not compr:
                data = fileread.read_gambit_neutral(fn)
            elif ftype in [ 'vtp', 'vtk'] and not compr:
                return read_vtk_surface(fn)
            else:
                raise "Unknown TriSurface type, cannot read file %s" % fn
            return TriSurface(*data)

        def write(self,fname,ftype=None,color=None):
            """Write the surface to file.

            If no filetype is given, it is deduced from the filename extension.
            If the filename has no extension, the 'off' file type is used.
            For a file with extension 'stl', the ftype may be 'stla' or 'stlb'
            to force ascii or binary STL format.
            The color is only useful for 'stlb' format.
            """
            if ftype is None:
                ftype,compr = utils.fileTypeComprFromExt(fname)
            else:
                ftype,compr = utils.fileTypeComprFromExt('a.'+ftype)

            if compr:
                raise ValueError("Compressed surface export is not active (yet)")

            print("Writing surface to file %s (%s)" % (fname, ftype))
            if ftype == 'pgf':
                Geometry.write(self, fname)
            elif ftype == 'gts':
                filewrite.writeGTS(fname, self.coords, self.getEdges(), self.getElemEdges())
                print("Wrote %s vertices, %s edges, %s faces" % self.shape())
            elif ftype in ['stl', 'stla', 'stlb', 'off', 'smesh', 'vtp', 'vtk']:
                if ftype in ['stl', 'stla']:
                    filewrite.writeSTL(fname, self.coords[self.elems], binary=False)
                elif ftype == 'stlb':
                    filewrite.writeSTL(fname, self.coords[self.elems], binary=True, color=color)
                elif ftype == 'off':
                    filewrite.writeOFF(fname, self)
                elif ftype == 'smesh':
                    from pyformex.plugins import tetgen
                    tetgen.writeSurface(fname, self.coords, self.elems)
                elif ftype == 'vtp' or ftype == 'vtk':
                    from pyformex.plugins import vtk_itf
                    vtk_itf.writeVTP(fname, self, checkMesh=False)
                print("Wrote %s vertices, %s elems" % (self.ncoords(), self.nelems()))
            else:
                print("Cannot save TriSurface as file %s" % fname)

    ####################### TriSurface Data ######################

        def avgVertexNormals(self):
            """Compute the average normals at the vertices."""
            return geomtools.averageNormals(self.coords, self.elems, atNodes=True)

        def areaNormals(self):
            """Compute the area and normal vectors of the surface triangles.

            The normal vectors are normalized.
            The area is always positive.

            The values are returned and saved in the object.
            """
            if self._areas is None or self._fnormals is None:
                self._areas, self._fnormals = geomtools.areaNormals(self.coords[self.elems])
            return self._areas, self._fnormals

        def areas(self):
            """Return the areas of all facets"""
            return self.areaNormals()[0]

        def volume(self):
            """Return the enclosed volume of the surface.

            This will only be correct if the surface is a closed manifold.
            """
            x = self.coords[self.elems]
            return inertia.surface_volume(x).sum()

        def volumeInertia(self,density=1.0):
            """Return the inertia properties of the enclosed volume of the surface.

            The surface should be a closed manifold and is supposed to be
            the border of a volume of constant density 1.

            Returns an :class:`inertia.Inertia` instance with attributes

            - `mass`: the total mass (float)
            - `ctr`: the center of mass: float (3,)
            - `tensor`: the inertia tensor in the central axes: shape (3,3)

            This will only be correct if the surface is a closed manifold.

            See :meth:`inertia` for the inertia of the surface.

            Example:

            >>> from pyformex.simple import sphere
            >>> S = sphere(8)
            >>> I = S.volumeInertia()
            >>> print(I.mass)  # doctest: +ELLIPSIS
            4.1526...
            >>> print(I.ctr)
            [ 0.  0.  0.]
            >>> print(I.tensor)
            [[ 1.65  0.   -0.  ]
             [ 0.    1.65 -0.  ]
             [-0.   -0.    1.65]]

            """
            x = self.coords[self.elems]
            V,C,I = inertia.surface_volume_inertia(x)
            I = inertia.Tensor(I)
            I = inertia.Inertia(I,mass=V,ctr=C)
            I.mass *= density
            I *= density
            return I

        def curvature(self,neighbours=1):
            """Return the curvature parameters at the nodes.

            This uses the nodes that are connected to the node via a shortest
            path of 'neighbours' edges.
            Eight values are returned: the Gaussian and mean curvature, the
            shape index, the curvedness, the principal curvatures and the
            principal directions.
            """
            curv = curvature(self.coords, self.elems, self.getEdges(), neighbours=neighbours)
            return curv

        def inertia(self,volume=False,density=1.0):
            """Return inertia related quantities of the surface.

            This computes the inertia properties of the centroids of the
            triangles, using the triangle area as a weight. The result is
            therefore different from self.coords.inertia() and usually better
            suited for the surface, especially if the triangle areas differ a lot.

            Returns a tuple with the center of gravity, the principal axes of
            inertia, the principal moments of inertia and the inertia tensor.

            See also :meth:`volumeInertia`.
            """
            if volume:
                return self.volumeInertia(density=density)
            else:
                I = self.centroids().inertia(mass=self.areas())
                I.mass *= density
                I *= density
                return I

        def surfaceType(self):
            """Check whether the TriSurface is a manifold and if it's closed."""
            ncon = self.nEdgeConnected()
            maxcon = ncon.max()
            mincon = ncon.min()
            manifold = maxcon == 2
            closed = mincon == 2
            return manifold, closed, mincon, maxcon

        def borderEdges(self):
            """Detect the border elements of TriSurface.

            The border elements are the edges having less than 2 connected elements.
            Returns True where edge is on the border.
            """
            return self.nEdgeConnected() <= 1

        def borderEdgeNrs(self):
            """Returns the numbers of the border edges."""
            return where(self.nEdgeConnected() <= 1)[0]

        def borderNodeNrs(self):
            """Detect the border nodes of TriSurface.

            The border nodes are the vertices belonging to the border edges.
            Returns a list of vertex numbers.
            """
            border = self.getEdges()[self.borderEdgeNrs()]
            return unique(border)

        ####### MANIFOLD #################

        def isManifold(self):
            """Check whether the TriSurface is a manifold.

            A surface is a manifold if a small sphere exists that cuts the surface
            to a surface that can continously be deformed to an open disk.
            """
            return self.surfaceType()[0]

        # TODO: We should probably add optional sorting (# connections) here
        # TODO: this could be moved into Mesh.nonManifoldEdges if it works
        #       for all level 2 Meshes
        def nonManifoldEdges(self):
            """Return the non-manifold edges.

            Non-manifold edges are edges having more than two triangles
            connected to them.

            Returns the indices of the non-manifold edges in a TriSurface.
            """
            return where(self.nEdgeConnected()>2)[0]

        def nonManifoldEdgesFaces(self):
            """Return the non-manifold edges and faces.

            Returns a tuple of:

            - the list of edges that connect 3 or more faces,
            - the list of faces connected to any of these edges.
            """
            conn = self.edgeConnections()
            ed = (conn!=-1).sum(axis=1)>2
            fa = unique(conn[ed])
            return arange(len(ed))[ed], fa[fa!=-1]

        def isClosedManifold(self):
            """Check whether the TriSurface is a closed manifold."""
            stype = self.surfaceType()
            return stype[0] and stype[1]

        def isConvexManifold(self):
            """Check whether the TriSurface is a convex manifold."""
            return self.isManifold() and self.edgeSignedAngles().min()>=0.

        def removeNonManifold(self):
            """Remove the non-manifold edges.

            Removes the non-manifold edges by iteratively applying two
            :meth:`removeDuplicate` and :meth:`collapseEdge` until no edge
            has more than two connected triangles.

            Returns the reduced surface.
            """
            S = self.removeDuplicate()
            non_manifold_edges = self.nonManifoldEdges()
            while non_manifold_edges.any():
                print("# nonmanifold edges: %s" % len(non_manifold_edges))
                maxcon = S.nEdgeConnected().max()
                wmax = where(S.nEdgeConnected()==maxcon)[0]
                S = S.collapseEdge(wmax[0])
                S = S.removeDuplicate()
                non_manifold_edges = S.nonManifoldEdges()
            return S

        ###### BORDER ######################

        def checkBorder(self):
            """Return the border of TriSurface.

            Returns a list of connectivity tables. Each table holds the
            subsequent line segments of one continuous contour of the border
            of the surface.
            """
            border = self.getEdges()[self.borderEdges()]
            if len(border) > 0:
                return connectedLineElems(border)
            else:
                return []

        def border(self,compact=True):
            """Return the border(s) of TriSurface.

            The complete border of the surface is returned as a list
            of plex-2 Meshes. Each Mesh constitutes a continuous part
            of the border. By default, the Meshes are compacted.
            Setting compact=False will return all Meshes with the full
            surface coordinate sets. This is usefull for filling the
            border and adding to the surface.
            """
            ML = [ Mesh(self.coords, e) for e in self.checkBorder() ]
            if compact:
                ML = [ M.compact() for M in ML ]
            return ML

        def fillBorder(self,method='radial',dir=None,compact=True):
            """Fill the border areas of a surface to make it closed.

            Returns a list of surfaces, each of which fills a singly connected
            part of the border of the input surface. Adding these surfaces to
            the original will create a closed surface. The surfaces will have
            property values set above those used in the parent surface.
            If the surface is already closed, an empty list is returned.

            There are three methods: 'radial', 'planar' and 'border',
            corresponding to the methods of the :func:`fillBorder` function.
            """
            if self.prop is None:
                mprop = 1
            else:
                mprop = self.prop.max()+1
            return [ fillBorder(b, method, dir).setProp(mprop+i) for i, b in enumerate(self.border(compact=compact)) ]

        def close(self,method='radial',dir=None):
            """This method needs documentation!!!!"""
            border = self.fillBorder(method, dir, compact=method=='radial')
            if method == 'radial':
                return self.concatenate([self]+border)
            else:
                elems = concatenate([ m.elems for m in [self]+border ], axis=0)
                if self.prop is None:
                    prop = zeros(shape=self.nelems(), dtype=Int)
                else:
                    prop = self.prop
                prop = concatenate( [prop] + [ m.prop for m in border ])
                return TriSurface(self.coords, elems, prop=prop)

        def edgeCosAngles(self,return_mask=False):
            """Return the cos of the angles over all edges.

            The surface should be a manifold (max. 2 elements per edge).
            Edges adjacent to only one element get cosangles = 1.0.
            If return_mask == True, a second return value is a boolean array
            with the edges that connect two faces.

            As a side effect, this method also sets the area, normals,
            elem_edges and edges attributes.
            """
            # get connections of edges to faces
            conn = self.getElemEdges().inverse()
            ## # BV: The following gives the same results, but does not
            ## #     guarantee the edges attribute to be set
            ## conn1 = self.edgeConnections()
            ## diff = (conn1 != conn).sum()
            ## if diff > 0:
            ##     print "edgeConnections",conn1
            ##     print "getElemEdges.inverse",conn
            # Bail out if some edge has more than two connected faces
            if conn.shape[1] > 2:
                raise RuntimeError("The TriSurface is not a manifold")
            # get normals on all faces
            n = self.areaNormals()[1]
            # Flag edges that connect two faces
            conn2 = (conn >= 0).sum(axis=-1) == 2
            # get adjacent facet normals for 2-connected edges
            n = n[conn[conn2]]
            # compute cosinus of angles over 2-connected edges
            cosa = dotpr(n[:, 0], n[:, 1])
            # Initialize cosangles to all 1. values
            cosangles = ones((conn.shape[0],))
            # Fill in the values for the 2-connected edges
            cosangles[conn2] = cosa
            # Clip to the -1...+1. range
            cosangles = cosangles.clip(min=-1., max=1.)
            # Return results
            if return_mask:
                return cosangles, conn2
            else:
                return cosangles

        def edgeAngles(self):
            """Return the angles over all edges (in degrees). It is the angle (0 to 180) between 2 face normals."""
            return arccosd(self.edgeCosAngles())

        def edgeSignedAngles(self,return_mask=False):
            """Return the signed angles over all edges (in degrees). It is the angle (-180 to 180) between 2 face normals.

            Positive/negative angles are associated to convexity/concavity at that edge.
            The border edges attached to one triangle have angle 0.
            NB: The sign of the angle is relevant if the surface has fixed normals. Should
            this check be done?
            """
            # get connections of edges to faces
            conn = self.getElemEdges().inverse()
            if conn.shape[1] > 2:
                raise RuntimeError("The TriSurface is not a manifold")
            # get normals on all faces
            n = self.areaNormals()[1]
            # Flag edges that connect two faces
            conn2 = (conn >= 0).sum(axis=-1) == 2
            # get adjacent facet normals for 2-connected edges
            n = n[conn[conn2]]
            edg = self.coords[self.getEdges()]
            edg = edg[:, 1]-edg[:, 0]
            ang = geomtools.rotationAngle(n[:, 0], n[:, 1],m=edg[conn2], angle_spec=DEG)
            # Initialize signed angles to all 0. values
            sangles = ones((conn.shape[0],))
            sangles[conn2] = ang
            # Clip to the -180...+180. range
            sangles = sangles.clip(min=-180., max=180.)
            # Return results
            if return_mask:
                return sangles, conn2
            else:
                return sangles

        def edgeLengths(self):
            """Returns the lengths of all edges

            Returns an array with the length of all the edges in the
            surface. As a side effect, this stores the connectivities
            of the edges to nodes and the elements to edges in the
            attributes edges, resp. elem_edges.
            """
            edg = self.coords[self.getEdges()]
            return length(edg[:, 1]-edg[:, 0])

        def _compute_data(self):
            """Compute data for all edges and faces."""
            if hasattr(self, 'edglen'):
                return
            self.areaNormals()
            edglen = self.edgeLengths()
            facedg = edglen[self.getElemEdges()]
            peri = facedg.sum(axis=-1)
            edgmin = facedg.min(axis=-1)
            edgmax = facedg.max(axis=-1)
            altmin = 2*self._areas / edgmax
            aspect = edgmax/altmin
            _qual_equi = sqrt(sqrt(3.)) / 6.
            qual = sqrt(self._areas) / peri / _qual_equi
            self.edglen, self.facedg, self.peri, self.edgmin, self.edgmax, self.altmin, self.aspect, self.qual = edglen, facedg, peri, edgmin, edgmax, altmin, aspect, qual

        def perimeters(self):
            """Compute the perimeters of all triangles."""
            self._compute_data()
            return self.peri

        def quality(self):
            """Compute a quality measure for the triangle schapes.

            The quality of a triangle is defined as the ratio of the square
            root of its surface area to its perimeter relative to this same
            ratio for an equilateral triangle with the same area.  The quality
            is then one for an equilateral triangle and tends to zero for a
            very stretched triangle.
            """
            self._compute_data()
            return self.qual

        def aspectRatio(self):
            """Return the apect ratio of the triangles of the surface.

            The aspect ratio of a triangle is the ratio of the longest edge
            over the smallest altitude of the triangle.

            Equilateral triangles have the smallest edge ratio (2 over square root 3).
            """
            self._compute_data()
            return self.aspect

        def smallestAltitude(self):
            """Return the smallest altitude of the triangles of the surface."""
            self._compute_data()
            return self.altmin

        def longestEdge(self):
            """Return the longest edge of the triangles of the surface."""
            self._compute_data()
            return self.edgmax

        def shortestEdge(self):
            """Return the shortest edge of the triangles of the surface."""
            self._compute_data()
            return self.edgmin

        def stats(self):
            """Return a text with full statistics."""
            bbox = self.bbox()
            manifold, closed, mincon, maxcon = self.surfaceType()
            self._compute_data()
            area = self.area()
            qual = self.quality()
            s = """
    Size: %s vertices, %s edges and %s faces
    Bounding box: min %s, max %s
    Minimal/maximal number of connected faces per edge: %s/%s
    Surface is%s a%s manifold
    Facet area: min %s; mean %s; max %s
    Edge length: min %s; mean %s; max %s
    Shortest altitude: %s; largest aspect ratio: %s
    Quality: %s .. %s
    """ % ( self.ncoords(), self.nedges(), self.nfaces(),
            bbox[0], bbox[1],
            mincon, maxcon,
            {True:'',False:' not'}[manifold], {True:' closed',False:''}[closed],
            self._areas.min(), self._areas.mean(), self._areas.max(),
            self.edglen.min(), self.edglen.mean(), self.edglen.max(),
            self.altmin.min(), self.aspect.max(),
            qual.min(), qual.max(),
            )
            if manifold:
                angles = self.edgeAngles()
                # getAngles is currently removed
                # vangles = self.getAngles()
                if closed:
                    volume = self.volume()

                s += """Angle between adjacent facets: min: %s; mean: %s; max: %s
    """ % ( angles.min(), angles.mean(), angles.max())

            s += "Total area: %s; " % area
            if manifold and closed:
                s += "Enclosed volume: %s" % volume
            else:
                s += "No volume (not a closed manifold!)"
            return s

        def distanceOfPoints(self,X,return_points=False):
            """Find the distances of points X to the TriSurface.

            The distance of a point is either:
            - the closest perpendicular distance to the facets;
            - the closest perpendicular distance to the edges;
            - the closest distance to the vertices.

            X is a (nX,3) shaped array of points.
            If return_points = True, a second value is returned: an array with
            the closest (foot)points matching X.
            """
            from pyformex.timer import Timer
            t = Timer()
            # distance from vertices
            ind, dist = geomtools.closest(X, self.coords, return_dist=True)
            if return_points:
                points = self.coords[ind]
            print("Vertex distance: %s seconds" % t.seconds(True))

            # distance from edges
            Ep = self.coords[self.getEdges()]
            res = geomtools.edgeDistance(X, Ep, return_points) # OKpid, OKdist, (OKpoints)
            okE, distE = res[:2]
            closer = distE < dist[okE]
            #print okE,closer
            if closer.size > 0:
                dist[okE[closer]] = distE[closer]
                if return_points:
                    points[okE[closer]] = res[2][closer]
            print("Edge distance: %s seconds" % t.seconds(True))
            #print dist

            # distance from faces
            Fp = self.coords[self.elems]
            res = geomtools.faceDistance(X, Fp, return_points) # OKpid, OKdist, (OKpoints)
            okF, distF = res[:2]
            closer = distF < dist[okF]
            #print okF,closer
            if closer.size > 0:
                dist[okF[closer]] = distF[closer]
                if return_points:
                    points[okF[closer]] = res[2][closer]
            print("Face distance: %s seconds" % t.seconds(True))
            #print dist

            if return_points:
                return dist, points
            else:
                return dist

        def degenerate(self):
            """Return a list of the degenerate faces according to area and normals.

            A face is degenerate if its surface is less or equal to zero or the
            normal has a nan.

            Returns the list of degenerate element numbers in a sorted array.
            """
            return geomtools.degenerate(*self.areaNormals())

        def removeDegenerate(self,compact=False):
            """Remove the degenerate elements from a TriSurface.

            Returns a TriSurface with all degenerate elements removed.
            By default, the coords set is unaltered and will still contain
            all points, even ones that are no longer connected to any element.
            To reduce the coordinate set, set the compact argument to True
            or use the :meth:`compact` method afterwards.
            """
            return self.cselect(self.degenerate(), compact=compact)

        def collapseEdge(self,edg):
            """Collapse an edge in a TriSurface.

            Collapsing an edge removes all the triangles connected to the edge
            and replaces the two vertices by a single one, placed at the center
            of the edge. Triangles only having one of the edge vertices, will
            all be connected to the new vertex.

            Parameters:

            - `edg`: the index of the edg to be removed. This is an index in the
              array of edges as returned by :meth:`getElemEdges`.

            Returns a TriSurface with the specified edge number removed.
            """
            # remove the elements connected to the collapsing edge
            invee = self.getElemEdges().inverse()
            els = invee[edg]
            els = els[els>=0]
            keep = complement(els,n=self.nelems())
            elems = self.elems[keep]
            prop = self.prop
            if prop is not None:
                prop = prop[keep]
            # replace the first node with the mean of the two
            node0,node1 = self.edges[edg]
            #print("Collapsing nodes %s and %s" % (node0,node1))
            elems[elems==node1] = node0
            coords = self.coords.copy()
            coords[node0] = 0.5 * ( coords[node0] + coords[node1] )
            return TriSurface(coords,elems,prop=prop).compact()

    ##################  Transform surface #############################
        # All transformations now return a new surface

        def offset(self,distance=1.):
            """Offset a surface with a certain distance.

            All the nodes of the surface are translated over a specified distance
            along their normal vector.
            """
            n = self.avgVertexNormals()
            coordsNew = self.coords + n*distance
            return TriSurface(coordsNew, self.getElems(), prop=self.prop)

        def dualMesh(self, method='median'):
            """Return the dual mesh of a compacted triangulated surface.

            It creates a new triangular mesh where all triangles with prop `p`
            represent the dual mesh region around the original surface node `p`.
            For more info, see http://users.led-inc.eu/~phk/mesh-dualmesh.html.

            - `method`: 'median' or 'voronoi'.

            Returns:

            - `method` = 'median': the Median dual mesh and the area of the region
              around each node. The sum of the node-based areas is equal to the
              original surface area.
            - `method` = 'voronoi': the Voronoi polyeders and a None.
            """
            if self.ncoords()!=self.compact().ncoords():
                raise ValueError("Expected a compacted surface")
            Q = self.convert('quad4', fuse=False)
            if method == 'voronoi':
                from pyformex.geomtools import triangleCircumCircle
                Q.coords[-self.nelems():] = triangleCircumCircle(self.coords[self.elems], bounding=False)[1]
            nconn = Q.nodeConnections()[arange(self.ncoords())]
            p = zeros(Q.nelems(), dtype=int)
            for i, conn in enumerate(nconn):
                p[conn[conn>-1]]=i
            Q = Q.setProp(p)
            if method == 'voronoi':
                return Q, None
            nodalAreas = asarray([Q.selectProp(i).area() for i in range(len(Q.propSet()))])
            return Q, nodalAreas

    ##################  Partitioning a surface  #############################

        def featureEdges(self,angle=60.):
            """Return the feature edges of the surface.

            Feature edges are edges that are prominent features of the geometry.
            They are either border edges or edges where the normals on the two
            adjacent triangles differ more than a given angle.
            The non feature edges then represent edges on a rather smooth surface.

            Parameters:

            - `angle`: The angle by which the normals on adjacent triangles
              should differ in order for the edge to be marked as a feature.

            Returns a boolean array with shape (nedg,) where the feature angles
            are marked with True.

            .. note::

               As a side effect, this also sets the `elem_edges` and `edges`
               attributes, which can be used to get the edge data with the same
               numbering as used in the returned mask. Thus, the following
               constructs a Mesh with the feature edges of a surface S::

                 p = S.featureEdges()
                 Mesh(S.coords,S.edges[p])
            """
            # Get the edge angles
            cosangles, conn2 = self.edgeCosAngles(return_mask=True)
            # initialize all edges as features
            feature = ones((self.edges.shape[0],), dtype=bool)
            # unmark edges with small angle
            feature[conn2] = cosangles[conn2] <= cosd(angle)
            return feature

        def partitionByAngle(self,angle=60.,sort='number'):
            """Partition the surface by splitting it at sharp edges.

            The surface is partitioned in parts in which all elements can be
            reach without ever crossing a sharp edge angle. More precisely,
            any two elements that can be connected by a line not crossing an
            edge between two elements having their normals differ more than
            angle (in degrees), will belong to the same part.

            The partitioning is returned as an integer array specifying the part
            number for eacht triangle.

            By default the parts are assigned property numbers in decreasing
            order of the number of triangles in the part. Setting the sort
            argument to 'area' will sort the parts according to decreasing
            area. Any other value will return the parts unsorted.

            Beware that the existence of degenerate elements may cause
            unexpected results. If unsure, use the :meth:`removeDegenerate`
            method first to remove those elements.
            """
            feat = self.featureEdges(angle=angle)
            p = self.maskedEdgeFrontWalk(mask=~feat, frontinc=0)

            if sort == 'number':
                p = sortSubsets(p)
            elif sort == 'area':
                p = sortSubsets(p, self.areaNormals()[0])

            return p

        # This may replace CutWithPlane after it has been proved stable
        # and has been expanded to multiple planes
        def cutWithPlane1(self,p,n,side='',return_intersection=False,atol=0.):
            """Cut a surface with a plane.

            Cuts the surface with a plane defined by a point p and normal n.

            Parameters:

            - `p`: float, shape (3,): a point in the cutting plane
            - `n`: float, shape (3,): the normal vector to the plane
            - `side`: '', '+' or '-': selector of the returned parts. Default
              is to return a tuple of two surfaces, with the parts at the positive,
              resp. negative side of the plane as defined by the normal vector.
              If a '+' or '-' is specified, only the corresponding part is returned.

            Returns:

            A tuple of two TriSurfaces, or a single TriSurface,
            depending on the value of `side`. The returned surfaces will have
            their normals fixed wherever possible. Property values will be set
            containing the triangle number of the original surface from which
            the elements resulted.
            """
            def finalize(Sp, Sn, I):
                # Result
                res = []
                if side in '+':
                    Sp = Sp.compact()#.fixNormals()
                    res.append(Sp)
                if side in '-':
                    Sn = Sn.compact()#.fixNormals()
                    res.append(Sn)
                if return_intersection:
                    res.append(I)
                if len(res) == 1:
                    res = res[0]
                else:
                    res = tuple(res)
                return res

            from pyformex.formex import _sane_side, _select_side
            side = _sane_side(side)

            try:
                p = array(p).reshape(3)
                n = array(n).reshape(3)
            except:
                raise ValueError("Expected a (3) shaped float array for both `p` and `n`")

            # Make sure we inherit element number
            save_prop = self.prop
            self.prop = arange(self.elems.shape[0])

            # Compute distance to plane of all vertices
            d = self.distanceFromPlane(p, n)

            p_pos = d > 0.
            p_neg = d < 0.
            p_in = ~(p_pos+p_neg)
            p_posin = p_pos + p_in
            p_negin = p_neg + p_in

            # Reduce the surface to the part intersecting with the plane:
            # Remember triangles with all vertices at same side
            # Elements completely in the plane end up in both parts.
            # BV: SHOULD WE CHANGE THIS???
            # TODO: put them only in negative?, as the volume is at the
            # negative side of the elements.
            all_p = p_posin[self.elems].all(axis=-1)
            all_n = p_negin[self.elems].all(axis=-1)
            S = self.cclip(all_p+all_n, compact=False)  # DOES THIS COMPACT? NO
            Sp = self.clip(all_p, compact=False)
            Sn = self.clip(all_n, compact=False)
            # Restore properties
            self.prop = save_prop

            #print "POS: %s; NEG: %s; CUT: %s" % (Sp.nelems(),Sn.nelems(),S.nelems())
            #clear()
            #drawPlane(p,n,((4.,4.),(4.,4.)))
            #draw(S,color='green')
            #drawNumbers(S.coords)

            # If there is no intersection, we're done
            # (we might have cut right along facet edges!)
            if S.nelems() == 0:
                res = _select_side(side, [Sp, Sn])
                return res

            #######################
            # Cut S with the plane.
            #######################
            # First define facets in terms of edges
            coords = S.coords
            edg = S.getEdges()
            fac = S.getElemEdges()
            ele = S.elems

            # Get the edges intersecting with the plane: 1 up and 1 down vertex
            d_edg = d[edg]
            edg_1_up = (d_edg > 0.).sum(axis=1) == 1
            edg_1_do = (d_edg < 0.).sum(axis=1) == 1
            cutedg = edg_1_up * edg_1_do
            ind = where(cutedg)[0]
            if ind.size == 0:
                raise ValueError("This really should not happen!")

            # Compute the intersection points
            M = Mesh(S.coords, edg[cutedg])
            x = geomtools.intersectionPointsSWP(M.toFormex().coords, p, n, mode='pair', return_all=True, atol=atol).reshape(-1, 3)
            # Create inverse index to lookup the point using the edge number
            rev = inverseUniqueIndex(ind) + coords.shape[0]
            # Concatenate the coords arrays
            coords = coords.concatenate([coords, x])

            # For each triangle, compute the number of cutting edges
            cut = cutedg[fac]
            ncut = cut.sum(axis=1)
            #print 'ncut',ncut

            if (ncut < 1).any() or (ncut > 2).any():
                # Maybe we should issue a warning and ignore these cases?
                print("NCUT: ", ncut)
                raise ValueError("I expected all triangles to be cut along 1 or 2 edges. I do not know how to proceed now.")

            if return_intersection:
                I = Mesh(eltype='line2')

            # Process the elements cutting one edge
            #######################################
            #print "Cutting 1 edge"
            ncut1 = ncut==1
            if ncut1.any():
                prop1 = where(ncut1)[0]
                fac1 = fac[ncut1]
                ele1 = ele[ncut1]
                #print ele1

                cutedg1 = cutedg[fac1]
                cut_edges =  fac1[cutedg1]

                # Identify the node numbers
                # 0 : vertex on positive side
                # 1 : vertex in plane
                # 2 : new point dividing edge
                # 3 : vertex on negative side
                elems = column_stack([
                    ele1[p_pos[ele1]],
                    ele1[p_in[ele1]],
                    rev[cut_edges],
                    ele1[p_neg[ele1]],
                    ])

                if side in '+':
                    Sp += TriSurface(coords, elems[:, 0:3], prop=prop1)
                if side in '-':
                    Sn += TriSurface(coords, elems[:, 1:4], prop=prop1)

            # Process the elements cutting two edges
            ########################################
            print("Cutting 2 edges")
            ncut2 = ncut==2     # selector over whole range
            print(ncut)
            print(ncut2)
            print(p_pos.sum(axis=-1)==2)
            if ncut2.any():
                prop2 = where(ncut2)[0]
                fac2 = fac[ncut2]
                ele2 = ele[ncut2]
                pp2 = p_pos[ele2]
                print("ele", ele2, pp2)
                ncut2p = pp2.sum(axis=-1)==1   # selector over ncut2 elems
                ncut2n = pp2.sum(axis=-1)==2
                print(ncut2p, ncut2n)

                if ncut2p.any():
                    #print "# one vertex at positive side"
                    prop1 = prop2[ncut2p]
                    fac1 = fac2[ncut2p]
                    ele1 = ele2[ncut2p]

                    print("ele1,1p", ele1)
                    cutedg1 = cutedg[fac1]
                    print(cutedg, fac1, cutedg1, fac1[cutedg1])
                    cut_edges =  fac1[cutedg1].reshape(-1, 2)
                    #print cut_edges

                    corner = ele1[p_pos[ele1]]
                    #print corner
                    quad = edg[cut_edges].reshape(-1, 4)
                    #print quad
                    #print quad != corner.reshape(-1,1)
                    quad2 = quad[ quad != corner.reshape(-1, 1) ]
                    #print quad2
                    # Identify the node numbers
                    # 0   : vertex on positive side
                    # 1,2 : new points dividing edges
                    # 3,4 : vertices on negative side
                    elems = column_stack([
                        ele1[p_pos[ele1]],
                        rev[cut_edges],
                        quad2.reshape(-1, 2)
                        # ele1[p_neg[ele1]].reshape(-1,2),
                        ])
                    #print elems

                    if side in '+':
                        Sp += TriSurface(coords, elems[:, 0:3], prop=prop1)
                    if side in '-':
                        Sn += TriSurface(coords, elems[:, 1:4], prop=prop1)
                        Sn += TriSurface(coords, elems[:, 2:5], prop=prop1)

                if ncut2n.any():
                    #print "# one vertex at negative side"
                    prop1 = where(ncut2n)[0]
                    fac1 = fac[ncut2n]
                    ele1 = ele[ncut2n]

                    cutedg1 = cutedg[fac1]
                    cut_edges =  fac1[cutedg1].reshape(-1, 2)
                    #print cut_edges

                    corner = ele1[p_neg[ele1]]
                    #print corner
                    quad = edg[cut_edges].reshape(-1, 4)
                    #print quad
                    #print quad != corner.reshape(-1,1)
                    quad2 = quad[ quad != corner.reshape(-1, 1) ]
                    #print quad2

                    # 0   : vertex on negative side
                    # 1,2 : new points dividing edges
                    # 3,4 : vertices on positive side
                    elems = column_stack([
                        quad2.reshape(-1, 2),
                        # we can not use this, because order of the 2 vertices
                        # is importtant
                        # ele1[p_pos[ele1]].reshape(-1,2),
                        rev[cut_edges],
                        ele1[p_neg[ele1]],
                        ])
                    #print elems

                    if side in '+':
                        Sp += TriSurface(coords, elems[:, 0:3], prop=prop1)
                        Sp += TriSurface(coords, elems[:, 1:4], prop=prop1)
                    if side in '-':
                        Sn += TriSurface(coords, elems[:, 2:5], prop=prop1)

            return finalize(Sp, Sn, I)
            # Result
            if side in '+':
                Sp = Sp.compact()#.fixNormals()
            if side in '-':
                Sn = Sn.compact()#.fixNormals()
            return _select_side(side, [Sp, Sn])

        def cutWithPlane(self,*args,**kargs):
            """Cut a surface with a plane or a set of planes.

            Cuts the surface with one or more plane and returns either one side
            or both.

            Parameters:

            - `p`,`n`: a point and normal vector defining the cutting plane.
              p and n can be sequences of points and vector,
              allowing to cut with multiple planes.
              Both p and n have shape (3) or (npoints,3).

            The parameters are the same as in :meth:`Formex.CutWithPlane`.
            The returned surface will have its normals fixed wherever possible.
            """
            F = self.toFormex()
            R = F.cutWithPlane(*args,**kargs)
            if isinstance(R, list):
                return [ TriSurface(r).fixNormals() for r in R ]
            else:
                return TriSurface(R).fixNormals()

    # TODO:
        ## UNDOCUMENTED! BECAUSE OF BAD FORMATTING
        ##
        ## LARGE CODE EXAMPLES SHOULD NOT GO IN THE DOCSTRINGS
        ## BUT IN AN EXAMPLE
        ##
        ## This should be in geomtools
        ## if it is better than  geomtools.intersectionPointsLWT,
        ## it sohuld replace oit.
    ##

        def intersectionWithLines(self,q,q2, method='line',atol=1.e-5):
            """Intersects a surface with lines.

            Parameters:

            - `q`,`q2`: (...,3) shaped arrays of points, defining
              a set of lines.
            - `method`: a string (line, segment or ray) defining if the line
              is either a full-line or a line-segment (q-q2) or a line-ray (q->q2)
            - `atol` : detected intersection points on the border edges (otherwise
              geomtools.insideTriangle could fail)

            Returns
            -------
            - a fused set of intersection points (Coords)
            - a (1,3) array with the indices of intersection point, line and
              triangle

            If a line is laying (parallel) on a triangle it will not generate
            intersections.
            There is a similar function in geomtools.intersectionPointsLWT but
            this one seems faster.
            """

            def insideLine(v,v0,v1,atol, method='segment'):
                """
                Check points inside a segments or a ray

                if method=='segment' check points inside segments v0,v1
                if method=='ray' check points inside ray v0->v1
                """
                if method=='segment':
                    ir=length(v-v0)+ length(v-v1)< length(v1-v0)+atol
                if method=='ray':
                    ir0 = dotpr(v-v0, normalize(v1-v0))>0#equivalent to ir0=length(normalize(v-v0)-m)<atol
                    ir1 = length(v-v0)<atol#point close to the ray's end
                    ir=ir0+ir1
                return ir

            r, C, n = geomtools.triangleBoundingCircle(self.coords[self.elems])#triangles' bounding sphere
            m = normalize(q2-q)
            ip = geomtools.pointNearLine(C, q, m, r, -1)#detects candidate lines/triangles (slow part)
            il = [ [ i ] * len(ip[i]) for i in range(len(ip)) ]
            t, l = concatenate(ip),concatenate(il).astype(Int)
            p = geomtools.intersectionPointsLWP(q[l], m[l], C[t], self.areaNormals()[1][t],  mode='pair')#intersects candidate lines/triangles
            prl=where(sum(isinf(p) + isnan(p), axis=1)>0)[0]#remove nan/inf (lines parallel to triangles)
            i1 = complement(prl, len(p))
            if len(i1)==0:
                return Coords(), []
            #
            # TODO
            # The use of shrink here is not a good technique!!!
            # Better is to include a tolerance in geomtools.insideTriangle
            #
            xt = self.select(t).toFormex().shrink(1.+atol)[:]#atol: insideTriangle sometimes fails on border!!
            xt, xp, xl = xt[i1], p[i1], l[i1]
            i2 = geomtools.insideTriangle(xt, xp[newaxis, ...]).reshape(-1)#remove intersections outside triangles
            i=i1[i2]
            if method!='line':
                xp, xl=xp[i2], xl[i2]
                i3=insideLine(xp, q[xl], q2[xl], atol, method)
                i=i[i3]
            p, j=p[i].fuse()
            return p, column_stack([j, l[i], t[i]])

        def intersectionWithPlane(self,p,n,atol=0.,sort='number'):
            """Return the intersection lines with plane (p,n).

            Returns a plex-2 mesh with the line segments obtained by cutting
            all triangles of the surface with the plane (p,n)
            p is a point specified by 3 coordinates.
            n is the normal vector to a plane, specified by 3 components.

            The return value is a plex-2 Mesh where the line segments defining
            the intersection are sorted to form continuous lines. The Mesh has
            property numbers such that all segments forming a single continuous
            part have the same property value.

            By default the parts are assigned property numbers in decreasing
            order of the number of line segments in the part. Setting the sort
            argument to 'distance' will sort the parts according to increasing
            distance from the point p.

            The splitProp() method can be used to get a list of Meshes.
            """
            n = asarray(n)
            p = asarray(p)
            # The vertices are classified based on their distance d from the plane:
            # - inside: d = 0
            # - up: d > 0
            # - down: d < 0

            # First, reduce the surface to the part intersecting with the plane:
            # remove triangles with all up or all down vertices
            d = self.distanceFromPlane(p, n)
            d_ele = d[self.elems]
            ele_all_up = (d_ele > 0.).all(axis=1)
            ele_all_do = (d_ele < 0.).all(axis=1)
            S = self.cselect(ele_all_up+ele_all_do, compact=False)

            # If there is no intersection, we're done
            if S.nelems() == 0:
                return Mesh(Coords(), Connectivity(nplex=2, eltype='line2'))

            Mparts = []
            coords = S.coords
            edg = S.getEdges()
            fac = S.getElemEdges()
            ele = S.elems
            # No need to recompute distances, as clipping does not compact!
            #d = S.distanceFromPlane(p,n) #

            # Get the edges intersecting with the plane: 1 up and 1 down vertex
            d_edg = d[edg]
            edg_1_up = (d_edg > 0.).sum(axis=1) == 1
            edg_1_do = (d_edg < 0.).sum(axis=1) == 1
            w = edg_1_up * edg_1_do
            ind = where(w)[0]

            # compute the intersection points
            if ind.size != 0:
                rev = inverseUniqueIndex(ind)
                M = Mesh(S.coords, edg[w])
                x = geomtools.intersectionPointsSWP(M.toFormex().coords, p, n, mode='pair', return_all=True, atol=atol).reshape(-1, 3)

            # For each triangle, compute the number of cutting edges
            cut = w[fac]
            ncut = cut.sum(axis=1)

            # Split the triangles based on the number of inside vertices
            d_ele = d[ele]
            ins = d_ele == 0.
            nins = ins.sum(axis=1)
            ins0, ins1, ins2, ins3 = [ where(nins==i)[0] for i in range(4) ]

            # No inside vertices -> 2 cutting edges
            if ins0.size > 0:
                cutedg = fac[ins0][cut[ins0]].reshape(-1, 2)
                e0 = rev[cutedg]
                Mparts.append(Mesh(x, e0, eltype='line2').compact())

            # One inside vertex
            if ins1.size > 0:
                ncut1 = ncut[ins1]
                cut10, cut11 = [ where(ncut1==i)[0] for i in range(2) ]
                # 0 cutting edges: does not generate a line segment
                # 1 cutting edge
                if cut11.size != 0:
                    e11ins = ele[ins1][cut11][ins[ins1][cut11]].reshape(-1, 1)
                    cutedg = fac[ins1][cut11][cut[ins1][cut11]].reshape(-1, 1)
                    e11cut = rev[cutedg]
                    x11 = Coords.concatenate([coords, x], axis=0)
                    e11 = concatenate([e11ins, e11cut+len(coords)], axis=1)
                    Mparts.append(Mesh(x11, e11, eltype='line2').compact())

            # Two inside vertices -> 0 cutting edges
            if ins2.size > 0:
                e2 = ele[ins2][ins[ins2]].reshape(-1, 2)
                Mparts.append(Mesh(coords, e2, eltype='line2').compact())

            # Three inside vertices -> 0 cutting edges
            if ins3.size > 0:
                insedg =  fac[ins3].reshape(-1)
                insedg.sort(axis=0)
                double = insedg == roll(insedg, 1, axis=0)
                insedg = setdiff1d(insedg, insedg[double])
                if insedg.size != 0:
                    e3 = edg[insedg]
                    Mparts.append(Mesh(coords, e3, eltype='line2').compact())

            # Done with getting the segments
            if len(Mparts) ==  0:
                # No intersection: return empty mesh
                return Mesh(Coords(), Connectivity(nplex=2, eltype='line2'))
            else:
                M = Mesh.concatenate(Mparts)

                # Remove degenerate and duplicate elements
                M = Mesh(M.coords, M.elems.removeDegenerate().removeDuplicate())

                # Split in connected loops
                parts = connectedLineElems(M.elems)
                prop = concatenate([ [i]*part.nelems() for i, part in enumerate(parts)])
                elems = concatenate(parts, axis=0)
                if sort == 'distance':
                    d = array([ M.coords[part].distanceFromPoint(p).min() for part in parts ])
                    srt = argsort(d)
                    inv = inverseUniqueIndex(srt)
                    prop = inv[prop]
                return Mesh(M.coords, elems, prop=prop)

        def slice(self,dir=0,nplanes=20):
            """Intersect a surface with a sequence of planes.

            A sequence of nplanes planes with normal dir is constructed
            at equal distances spread over the bbox of the surface.

            The return value is a list of intersectionWithPlane() return
            values, i.e. a list of Meshes, one for every cutting plane.
            In each Mesh the simply connected parts are identified by
            property number.
            """
            o = self.center()
            if isinstance(dir, int):
                dir = unitVector(dir)
            xmin, xmax = self.coords.directionalExtremes(dir, o)
            P = Coords.interpolate(xmin, xmax, nplanes)
            return [ self.intersectionWithPlane(p, dir) for p in P ]

    ##################  Smooth a surface #############################

        def smooth(self,method='lowpass',iterations=1,lambda_value=0.5,neighbourhood=1,alpha=0.0,beta=0.2):
            """Smooth the surface.

            Returns a TriSurface which is a smoothed version of the original.
            Two smoothing methods are available: 'lowpass' and 'laplace'.

            Parameters:

            - `method`: 'lowpass' or 'laplace'
            - `iterations`: int: number of iterations
            - `lambda_value`: float: lambda value used in the filters

            Extra parameters for 'lowpass' and 'laplace':

            - `neighbourhood`: int: maximum number of edges followed in defining
              the node neighbourhood

            Extra parameters for 'laplace':

            - `alpha`, `beta`: float: parameters for the laplace method.

            Returns the smoothed TriSurface
            """
            method = method.lower()

            # find adjacency
            adj = adjacencyArrays(self.getEdges(), nsteps=neighbourhood)[1:]
            adj = column_stack(adj)
            # find interior vertices
            bound_edges = self.borderEdgeNrs()
            inter_vertex = resize(True, self.ncoords())
            inter_vertex[unique(self.getEdges()[bound_edges])] = False
            # calculate weights
            w = ones(adj.shape, dtype=float)
            w[adj<0] = 0.
            val = (adj>=0).sum(-1).reshape(-1, 1)
            w /= val
            w = w.reshape(adj.shape[0], adj.shape[1], 1)

            # recalculate vertices

            if method == 'laplace':
                xo = self.coords
                x = self.coords.copy()
                for step in range(iterations):
                    xn = x + lambda_value*(w*(x[adj]-x.reshape(-1, 1, 3))).sum(1)
                    xd = xn - (alpha*xo + (1-alpha)*x)
                    x[inter_vertex] = xn[inter_vertex] - (beta*xd[inter_vertex] + (1-beta)*(w[inter_vertex]*xd[adj[inter_vertex]]).sum(1))

            else: # default: lowpass
                k = 0.1
                mu_value = -lambda_value/(1-k*lambda_value)
                x = self.coords.copy()
                for step in range(iterations):
                    x[inter_vertex] = x[inter_vertex] + lambda_value*(w[inter_vertex]*(x[adj[inter_vertex]]-x[inter_vertex].reshape(-1, 1, 3))).sum(1)
                    x[inter_vertex] = x[inter_vertex] + mu_value*(w[inter_vertex]*(x[adj[inter_vertex]]-x[inter_vertex].reshape(-1, 1, 3))).sum(1)

            return TriSurface(x, self.elems, prop=self.prop)

        def smoothLowPass(self,iterations=2,lambda_value=0.5,neighbours=1):
            """Apply a low pass smoothing to the surface."""
            return self.smooth('lowpass', iterations//2, lambda_value, neighbours)

        def smoothLaplaceHC(self,iterations=2,lambda_value=0.5,alpha=0.,beta=0.2):
            """Apply Laplace smoothing with shrinkage compensation to the surface."""
            return self.smooth('laplace', iterations, lambda_value, alpha=alpha, beta=beta)

        def refine(self,max_edges=None,min_cost=None,method='gts'):
            """Refine the TriSurface.

            Refining a TriSurface means increasing the number of triangles and
            reducing their size, while keeping the changes to the modeled surface
            minimal.
            Construct a refined version of the surface.
            This uses the external program `gtsrefine`. The surface
            should be a closed orientable non-intersecting manifold.
            Use the :meth:`check` method to find out.

            Parameters:

            - `max_edges`: int: stop the refining process if the number of
              edges exceeds this value
            - `min_cost`: float: stop the refining process if the cost of refining
              an edge is smaller
            - `log`: boolean: log the evolution of the cost
            - `verbose`: boolean: print statistics about the surface
            """
            if method == 'gts':
                return self.gts_refine(max_edges, min_cost)

            # THIS IS WORK IN PROGRESS
            self.getElemEdges()
            edglen = length(self.coords[self.edges[:, 1]]-self.coords[self.edges[:, 0]])
            print(edglen)
            return self

        def similarity(self,S):
            """Compute the similarity with another TriSurface.

            Compute a quantitative measure of the similarity of the volumes
            enclosed by two TriSurfaces. Both the calling and the passed
            TriSurface should be closed manifolds (see :meth:`isClosedManifold`).

            Returns a tuple a tuple (jaccard, dice, overlap).
            If A and B are two closed manifolds, VA and VB are their respective
            volumes, VC is the volume of the intersection of A and B, and VD is
            the volume of the union of A and B, then the following similarity
            measures are defined:

            - jaccard coefficient: VC / VD
            - dice: 2 * VC / (VA + VB)
            - overlap: VC / min(VA,VB)

            Both jaccard and dice range from 0 when the surfaces are completely
            disjoint to 1 when the surfaces are identical. The overlap coefficient
            becomes 1 when one of the surfaces is completely inside the other.

            This method uses gts library to compute the intersection or union.
            If that fails, nan values are returned.
            """
            A,B = self,S
            VA = A.volume()
            VB = B.volume()
            try:
                VC = A.boolean(B,'*').volume()
                VD = VA + VB - VC
            except:
                try:
                    VD = A.boolean(B,'+').volume()
                    VC = VA + VB - VD
                except:
                    VC = VD = nan

            dice = 2 * VC / (VA+VB)
            overlap = VC / min([VA,VB])
            jaccard = VC / VD
            return  jaccard, dice, overlap

    ###################################################################
    ##    Methods using admesh/GTS
    ##############################

        def fixNormals(self,outwards=True):
            """Fix the orientation of the normals.

            Some surface operations may result in improperly oriented normals,
            switching directions from one triangle to the adjacent one.
            This method tries to reverse improperly oriented normals so that a
            singly oriented surface is achieved.

            .. note: In the current version, this uses the external program
               `admesh`, so this should be installed on the machine.

            If the surface is a (possibly non-orientable) manifold, the result
            will be an orientable manifold.

            If the surface is a closed manifold, the normals will be
            oriented to the outside. This is done by computing the volume
            inside the surface and reversing the normals if that turns out
            to be negative.

            Parameters:

            - `outwards`: boolean: if True (default), a test is done whether
              the surface is a closed manifold, and if so, the normals are
              oriented outwards. Setting this value to False will skip this
              test and the (possible) reversal of the normals.
            """
            if self.nelems() == 0:
                # Protect against impossible file handling for empty surfaces
                return self
            tmp = utils.tempName('.stl')
            tmp1 = utils.tempName('.off')
            print("Writing temp file %s" % tmp)
            self.write(tmp, 'stl')
            print("Fixing surface while converting to OFF format %s" % tmp1)
            stlConvert(tmp, tmp1)
            print("Reading result from %s" % tmp1)
            S = TriSurface.read(tmp1)
            os.remove(tmp)
            os.remove(tmp1)
            S.setProp(self.prop)

            if outwards and S.isClosedManifold() and S.volume() < 0.:
                S = S.reverse()

            return S

        def check(self,matched=True,verbose=False):
            """Check the surface using gtscheck.

            Uses `gtscheck` to check whether the surface is an orientable,
            non self-intersecting manifold.

            This is a necessary condition for using the `gts` methods:
            split, coarsen, refine, boolean. (Additionally, the surface should be
            closed, wich can be checked with :meth:`isClosedManifold`).

            Returns a tuple of:

            - an integer return code with the value:

              - 0: the surface is an orientable, non self-intersecting manifold.
              - 1: the created GTS file is invalid: this should normally not occur.
              - 2: the surface is not an orientable manifold. This may be due to
                misoriented normals. The :meth:`fixNormals` and :meth:`reverse`
                methods may be used to help fixing the problem in such case.
              - 3: the surface is an orientable manifold but is
                self-intersecting. The self intersecting triangles are returned as
                the second return value.

            - the intersecting triangles in the case of a return code 3, else None.
              If matched==True, intersecting triangles are returned as element
              indices of self, otherwise as a separate TriSurface object.

            If verbose is True, prints the statistics reported by the gtscheck
            command.
            """
            tmp = utils.tempName('.gts')
            self.write(tmp, 'gts')
            P = utils.system("gtscheck -v",stdin=tmp)
            if verbose:
                print(P.sta)
            os.remove(tmp)
            if P.sta == 0:
                print('The surface is an orientable non self-intersecting manifold')
                return P.sta, None
            if P.sta==2:
                print('The surface is not an orientable manifold (this may be due to badly oriented normals)')
                return P.sta, None
            if P.sta==3:
                print('The surface is an orientable manifold but is self-intersecting')
                tmp = utils.tempName('.gts')
                print("Writing temp file %s" % tmp)
                fil = open(tmp, 'w')
                fil.write(P.out)
                fil.close()
                Si = TriSurface.read(tmp)
                os.remove(tmp)
                if matched:
                    return P.sta, self.matchCentroids(Si)
                else:
                    return P.sta, Si
            else:
                print('Status of gtscheck not understood')
                return P.sta, None

        def split(self,base,verbose=False):
            """Split the surface using gtssplit.

            Splits the surface into connected and manifold components.
            This uses the external program `gtssplit`. The surface
            should be a closed orientable non-intersecting manifold.
            Use the :meth:`check` method to find out.

            This method creates a series of files with given base name,
            each file contains a single connected manifold.
            """
            cmd = 'gtssplit -v %s' % base
            if verbose:
                cmd += ' -v'
            tmp = utils.tempName('.gts')
            print("Writing temp file %s" % tmp)
            self.write(tmp, 'gts')
            print("Splitting with command\n %s" % cmd)
            cmd += ' < %s' % tmp
            P = utils.command(cmd, shell=True)
            os.remove(tmp)
            if P.sta or verbose:
                print(P.out)
            #
            # TODO: WE SHOULD READ THIS FILES BACK !!!
            #

        def coarsen(self,min_edges=None,max_cost=None,
                    mid_vertex=False, length_cost=False, max_fold=1.0,
                    volume_weight=0.5, boundary_weight=0.5, shape_weight=0.0,
                    progressive=False, log=False, verbose=False):
            """Coarsen the surface using gtscoarsen.

            Construct a coarsened version of the surface.
            This uses the external program `gtscoarsen`. The surface
            should be a closed orientable non-intersecting manifold.
            Use the :meth:`check` method to find out.

            Parameters:

            - `min_edges`: int: stops the coarsening process if the number of
              edges was to fall below it
            - `max_cost`: float: stops the coarsening process if the cost of
              collapsing an edge is larger
            - `mid_vertex`: boolean: use midvertex as replacement vertex instead
              of the default, which is a volume optimized point
            - `length_cost`: boolean: use length^2 as cost function instead of the
              default optimized point cost
            - `max_fold`: float: maximum fold angle in degrees
            - `volume_weight`: float: weight used for volume optimization
            - `boundary_weight`: float: weight used for boundary optimization
            - `shape_weight`: float: weight used for shape optimization
            - `progressive`: boolean: write progressive surface file
            - `log`: boolean: log the evolution of the cost
            - `verbose`: boolean: print statistics about the surface
            """
            if min_edges is None and max_cost is None:
                min_edges = self.nedges() // 2
            cmd = 'gtscoarsen'
            if min_edges:
                cmd += ' -n %d' % min_edges
            if max_cost:
                cmd += ' -c %f' % max_cost
            if mid_vertex:
                cmd += ' -m'
            if length_cost:
                cmd += ' -l'
            if max_fold:
                cmd += ' -f %f' % max_fold
            cmd += ' -w %f' % volume_weight
            cmd += ' -b %f' % boundary_weight
            cmd += ' -s %f' % shape_weight
            if progressive:
                cmd += ' -p'
            if log:
                cmd += ' -L'
            if verbose:
                cmd += ' -v'
            tmp = utils.tempName('.gts')
            tmp1 = utils.tempName('.gts')
            print("Writing temp file %s" % tmp)
            self.write(tmp, 'gts')
            print("Coarsening with command\n %s" % cmd)
            P = utils.command(cmd, stdin=tmp, stdout=tmp1)
            os.remove(tmp)
            if P.sta or verbose:
                print(P.out)
            print("Reading coarsened model from %s" % tmp1)
            S = TriSurface.read(tmp1)
            os.remove(tmp1)
            return S

        def gts_refine(self,max_edges=None,min_cost=None,log=False,verbose=False):
            """Refine the TriSurface.

            Refining a TriSurface means increasing the number of triangles and
            reducing their size, while keeping the changes to the modeled surface
            minimal.
            Construct a refined version of the surface.
            This uses the external program `gtsrefine`. The surface
            should be a closed orientable non-intersecting manifold.
            Use the :meth:`check` method to find out.

            Parameters:

            - `max_edges`: int: stop the refining process if the number of
              edges exceeds this value
            - `min_cost`: float: stop the refining process if the cost of refining
              an edge is smaller
            - `log`: boolean: log the evolution of the cost
            - `verbose`: boolean: print statistics about the surface
            """
            if max_edges is None and min_cost is None:
                max_edges = self.nedges() * 2
            cmd = 'gtsrefine'
            if max_edges:
                cmd += ' -n %d' % max_edges
            if min_cost:
                cmd += ' -c %f' % min_cost
            if log:
                cmd += ' -L'
            if verbose:
                cmd += ' -v'
            tmp = utils.tempName('.gts')
            tmp1 = utils.tempName('.gts')
            print("Writing temp file %s" % tmp)
            self.write(tmp, 'gts')
            print("Refining with command\n %s" % cmd)
            P = utils.command(cmd, stdin=tmp, stdout=tmp1)
            os.remove(tmp)
            if P.sta or verbose:
                print(P.out)
            print("Reading refined model from %s" % tmp1)
            S = TriSurface.read(tmp1)
            os.remove(tmp1)
            return S

        def gts_smooth(self,iterations=1,lambda_value=0.5,verbose=False):
            """Smooth the surface using gtssmooth.

            Smooth a surface by applying iterations of a Laplacian filter.
            This uses the external program `gtssmooth`. The surface
            should be a closed orientable non-intersecting manifold.
            Use the :meth:`check` method to find out.

            Parameters:

            - `lambda_value`: float: Laplacian filter parameter
            - `iterations`: int: number of iterations
            - `verbose`: boolean: print statistics about the surface

            See also: :meth:`smoothLowPass`, :meth:`smoothLaplaceHC`
            """
            cmd = 'gtssmooth'
    #        if fold_smoothing:
    #            cmd += ' -f %s' % fold_smoothing
            cmd += ' %s %s' % (lambda_value, iterations)
            if verbose:
                cmd += ' -v'
            tmp = utils.tempName('.gts')
            tmp1 = utils.tempName('.gts')
            print("Writing temp file %s" % tmp)
            self.write(tmp, 'gts')
            print("Smoothing with command\n %s" % cmd)
            P = utils.command(cmd, stdin=tmp, stdout=tmp1)
            os.remove(tmp)
            if P.sta or verbose:
                print(P.out)
            print("Reading smoothed model from %s" % tmp1)
            S = TriSurface.read(tmp1)
            os.remove(tmp1)
            return S

        def inside(self,pts,method='gts',tol='auto',multi=False):
            """Test which of the points pts are inside the surface.

            Parameters:

            - `pts`: a Coords or compatible.
            - `method`: string: method to be used for the detection. Depending on
              the software you have installed the following are possible:

              - 'gts': provided by pyformex-extra
              - 'vtk': provided by python-vtk (slower)

            - `tol`: tolerance on equality of floating point values

            Returns an integer array with the indices of the points that are
            inside the surface. The indices refer to the onedimensional list
            of points as obtained from pts.points().
            """
            pts = Coords(pts).points()
            if method == 'gts':
                from pyformex.plugins.pyformex_gts import inside
                return inside(self, pts, tol, multi=multi)
            elif method == 'vtk':
                from pyformex.plugins.vtk_itf import inside
                return inside(self, pts, tol)

        def outside(self, pts, **kargs):
            return complement(self.inside(pts, **kargs), len(pts))

        def voxelize(self,n,bbox=0.01,return_formex=False):
            """Voxelize the volume inside a closed surface.

            Parameters:

            - `n`: int or (int, int, int): resolution, i.e. number of voxel cells
              to use along the three axes.
              If a single int is specified, the number of cells will be adapted
              according to the surface's :meth:`sizes` (as the voxel cells are
              always cubes). The specified number of voxels will be use along the
              largest direction.
            - `bbox`: float or (point,point): defines the bounding box of the volume
              that needs to be voxelized. A float specifies a relative amount to add
              to the surface's bounding box. Note that this defines the bounding box
              of the centers of the voxels.
            - `return_formex`: bool; if True, also returns a Formex with the centers
              of the voxels.

            Returns an int array of shape (nz,ny,nx) with value 1 for the voxels
            whose center is inside the surface, else 0.
            If `return_formex` is True, also returns a plex-1 Formex with the
            centers of the voxels , and property values 0 or 1 if the point is
            respectively outside or inside the surface. The voxel cell ordering
            in the Formex is z-direction first, then y, then x.

            See also example Voxelize, for saving the voxel values in a stack
            of binary images.

            """
            if not self.isClosedManifold():
                raise ValueError("The surface is non a closed manifold")

            from pyformex import simple
            if isFloat(bbox):
                a,b = 1.0+bbox, bbox
                bbox = self.bbox()
                bbox = [ a*bbox[0]-b*bbox[1], a*bbox[1]-b*bbox[0]]
            bbox = checkArray(bbox,shape=(2,3),kind='f')
            if isInt(n):
                sz = bbox[1]-bbox[0]
                step = sz.max() / (n-1)
                n = ceil(sz / step).astype(Int)
            n = checkArray(n,shape=(3,),kind='i')
            X = simple.regularGrid(bbox[0], bbox[0]+n*step, n, swapaxes=True)
            ind = self.inside(X)
            vox = zeros(n+1, dtype=uint8)
            vox.ravel()[ind] = 1
            if return_formex:
                P = Formex(X.reshape(-1,3))
                P.setProp(vox.ravel())
                return vox,P
            return vox

        def tetgen(self,quality=True,volume=None,filename=None,format='.off'):
            """Create a tetrahedral mesh inside the surface

            - `quality`: if True, the output will be a quality mesh
              The circumradius-to-shortest-edge ratio can be constrained by
              specifying a float value for quality (default is 2.0)
              - `volume`: float: applies a maximum tetrahedron volume constraint
            - `filename`: if specified, the surface model will be saved on this
              file and the tetgen models will be named likewise. If unspecified,
              temporary file names will be used.

            If the creation of the tetrahedral model is succesful, the
            resulting tetrahedral mesh is returned.
            """
            from pyformex.plugins import tetgen
            if filename is None:
                outputdir = utils.tempDir()
                fn = os.path.join(outputdir, 'surface') + format
            else:
                fn = filename
            self.write(fn)
            res = tetgen.tetMesh(fn, quality, volume)
            if filename is None:
                utils.removeTree(outputdir)
            return res['tetgen.ele']

        @utils.deprecated_by('TriSurface.facetArea','TriSurface.areas')
        def facetArea(self):
            return self.areas()

    ##########################################################################
    ################# Non-member and obsolete functions ######################

    def find_row(mat,row,nmatch=None):
        """Find all rows in matrix matching given row."""
        if nmatch is None:
            nmatch = mat.shape[1]
        return where((mat == row).sum(axis=1) == nmatch)[0]

    def find_nodes(nodes, coords):
        """Find nodes with given coordinates in a node set.

        nodes is a (nnodes,3) float array of coordinates.
        coords is a (npts,3) float array of coordinates.

        Returns a (n,) integer array with ALL the node numbers matching EXACTLY
        ALL the coordinates of ANY of the given points.
        """
        return concatenate([ find_row(nodes, c) for c in coords])

    def find_first_nodes(nodes, coords):
        """Find nodes with given coordinates in a node set.

        nodes is a (nnodes,3) float array of coordinates.
        coords is a (npts,3) float array of coordinates.

        Returns a (n,) integer array with THE FIRST node number matching EXACTLY
        ALL the coordinates of EACH of the given points.
        """
        res = [ find_row(nodes, c) for c in coords ]
        return array([ r[0] for r in res ])

    def find_triangles(elems, triangles):
        """Find triangles with given node numbers in a surface mesh.

        elems is a (nelems,3) integer array of triangles.
        triangles is a (ntri,3) integer array of triangles to find.

        Returns a (ntri,) integer array with the triangles numbers.
        """
        magic = elems.max()+1

        mag1 = enmagic3(elems, magic)
        mag2 = enmagic3(triangles, magic)

        nelems = elems.shape[0]
        srt = mag1.argsort()
        old = arange(nelems)[srt]
        mag1 = mag1[srt]
        pos = mag1.searchsorted(mag2)
        tri = where(mag1[pos]==mag2, old[pos], -1)
        return tri

    def remove_triangles(elems, remove):
        """Remove triangles from a surface mesh.

        elems is a (nelems,3) integer array of triangles.
        remove is a (nremove,3) integer array of triangles to remove.

        Returns a (nelems-nremove,3) integer array with the triangles of
        nelems where the triangles of remove have been removed.
        """
        print("Removing %s out of %s triangles" % (remove.shape[0], elems.shape[0]))
        magic = elems.max()+1

        mag1 = enmagic3(elems, magic)
        mag2 = enmagic3(remove, magic)

        mag1.sort()
        mag2.sort()

        nelems = mag1.shape[0]

        pos = mag1.searchsorted(mag2)
        mag1[pos] = -1
        mag1 = mag1[mag1 >= 0]

        elems = demagic3(mag1, magic)
        print("Actually removed %s triangles, leaving %s" % (nelems-mag1.shape[0], elems.shape[0]))

        return elems

    ####### Deprecated function #######################

    @utils.deprecated_by('trisurface.Sphere','simple.sphere')
    def Sphere(level=4,verbose=False,filename=None):
        from pyformex import simple
        return simple.sphere(ndiv=2**(level-1))

    @utils.deprecated_by('trisurface.intersectSurfaceWithLines','trisurface.intersectionWithLines')
    def intersectSurfaceWithLines(ts, qli, mli, atol=1.e-5):
        pass
    ##this function has been removed. You can get the same results using:
    ##      P,X=intersectionWithLines(self, q, m, q2,atol=1.e-5)
    ##      P,L,T = P[X[:, 0]], X[:, 1], X[:, 2]

    @utils.deprecated_by('trisurface.intersectSurfaceWithSegments','trisurface.intersectionWithLines')
    def intersectSurfaceWithSegments(s1, segm, atol=1.e-5):
        pass
    ##this function has been removed. You can get the same results using:
    ##      P,X=intersectionWithLines(self, q=segm[:, 0], m=None, q2=segm[:, 1],atol=1.e-5)
    ##      P,L,T = P[X[:, 0]], X[:, 1], X[:, 2]

    # Set TriSurface methods defined elsewhere

    from pyformex.plugins import pyformex_gts
    TriSurface.boolean = pyformex_gts.boolean
    TriSurface.intersection = pyformex_gts.intersection
    TriSurface.gtsset = pyformex_gts.gtsset

    from pyformex.plugins import webgl
    TriSurface.webgl = webgl.surface2webgl

    def read_vtk_surface(fn):
        try:
            from pyformex.plugins import vtk_itf
            return vtk_itf.read_vtk_surface(fn)
        except:
            utils.warn("I could not import VTK. This probably means that Python"
                       "VTK bindings are not installed on your machine.")

    # End