def smooth_mesh(mesh, n_iter=4, lam=0.6307, mu=-0.6347,
                weights=None, bconstr=True,
                volume_corr=False):
    """
    FE mesh smoothing.

    Based on:

    [1] Steven K. Boyd, Ralph Muller, Smooth surface meshing for automated
    finite element model generation from 3D image data, Journal of
    Biomechanics, Volume 39, Issue 7, 2006, Pages 1287-1295,
    ISSN 0021-9290, 10.1016/j.jbiomech.2005.03.006.
    (http://www.sciencedirect.com/science/article/pii/S0021929005001442)

    Parameters
    ----------
    mesh : mesh
        FE mesh.
    n_iter : integer, optional
        Number of iteration steps.
    lam : float, optional
        Smoothing factor, see [1].
    mu : float, optional
        Unshrinking factor, see [1].
    weights : array, optional
        Edge weights, see [1].
    bconstr: logical, optional
        Boundary constraints, if True only surface smoothing performed.
    volume_corr: logical, optional
        Correct volume after smoothing process.

    Returns
    -------
    coors : array
        Coordinates of mesh nodes.
    """

    def laplacian(coors, weights):

        n_nod = coors.shape[0]
        displ = (weights - sps.identity(n_nod)) * coors

        return displ

    def taubin(coors0, weights, lam, mu, n_iter):

        coors = coors0.copy()

        for ii in range(n_iter):
            displ = laplacian(coors, weights)
            if nm.mod(ii, 2) == 0:
                coors += lam * displ
            else:
                coors += mu * displ

        return coors

    def dets_fast(a):
        m = a.shape[0]
        n = a.shape[1]
        lapack_routine = lapack_lite.dgetrf
        pivots = nm.zeros((m, n), intc)
        flags = nm.arange(1, n + 1).reshape(1, -1)
        for i in xrange(m):
            tmp = a[i]
            lapack_routine(n, n, tmp, n, pivots[i], 0)
        sign = 1. - 2. * (nm.add.reduce(pivots != flags, axis=1) % 2)
        idx = nm.arange(n)
        d = a[:, idx, idx]
        absd = nm.absolute(d)
        sign *= nm.multiply.reduce(d / absd, axis=1)
        nm.log(absd, absd)
        logdet = nm.add.reduce(absd, axis=-1)

        return sign * nm.exp(logdet)

    def get_volume(el, nd):

        dim = nd.shape[1]
        nnd = el.shape[1]

        etype = '%d_%d' % (dim, nnd)
        if etype == '2_4' or etype == '3_8':
            el = elems_q2t(el)

        nel = el.shape[0]

        mul = 1.0 / factorial(dim)
        if dim == 3:
            mul *= -1.0

        mtx = nm.ones((nel, dim + 1, dim + 1), dtype=nm.double)
        mtx[:,:,:-1] = nd[el,:]
        vols = mul * dets_fast(mtx.copy()) # copy() ???
        vol = vols.sum()
        bc = nm.dot(vols, mtx.sum(1)[:,:-1] / nnd)

        bc /= vol

        return vol, bc


    n_nod = mesh.n_nod

    if weights is None:
        # initiate all vertices as inner - hierarchy = 2
        node_group = nm.ones((n_nod,), dtype=nm.int8) * 2
        sndi, edges = get_snodes_uedges(mesh.conns[0], mesh.descs[0])
        # boundary vertices - set hierarchy = 4
        if bconstr:
            node_group[sndi] = 4

        # generate costs matrix
        end1 = edges[:,0]
        end2 = edges[:,1]
        idxs = nm.where(node_group[end2] >= node_group[end1])
        rows1 = end1[idxs]
        cols1 = end2[idxs]
        idxs = nm.where(node_group[end1] >= node_group[end2])
        rows2 = end2[idxs]
        cols2 = end1[idxs]
        crows = nm.concatenate((rows1, rows2))
        ccols = nm.concatenate((cols1, cols2))
        costs = sps.coo_matrix((nm.ones_like(crows), (crows, ccols)),
                               shape=(n_nod, n_nod),
                               dtype=nm.double)

        # generate weights matrix
        idxs = range(n_nod)
        aux = sps.coo_matrix((1.0 / nm.asarray(costs.sum(1)).squeeze(),
                              (idxs, idxs)),
                             shape=(n_nod, n_nod),
                             dtype=nm.double)

        #aux.setdiag(1.0 / costs.sum(1))
        weights = (aux.tocsc() * costs.tocsc()).tocsr()

    coors = taubin(mesh.coors, weights, lam, mu, n_iter)

    if volume_corr:
        volume0, bc = get_volume(mesh.conns[0], mesh.coors)
        volume, _ = get_volume(mesh.conns[0], coors)

        scale = volume0 / volume
        coors = (coors - bc) * scale + bc

    return coors
def smooth_mesh ( mesh , n_iter = 4 , lam = 0.6307 , mu = - 0.6347 , weights = None , bconstr = True , volume_corr = False ) : def laplacian ( coors , weights ) : n_nod = coors . shape [ 0 ] displ = ( weights - sps . identity ( n_nod ) ) * coors return displ def taubin ( coors0 , weights , lam , mu , n_iter ) : coors = coors0 . copy ( ) for ii in range ( n_iter ) : displ = laplacian ( coors , weights ) if nm . mod ( ii , 2 ) == 0 : coors += lam * displ else : coors += mu * displ return coors def dets_fast ( a ) : m = a . shape [ 0 ] n = a . shape [ 1 ] lapack_routine = lapack_lite . dgetrf pivots = nm . zeros ( ( m , n ) , intc ) flags = nm . arange ( 1 , n + 1 ) . reshape ( 1 , - 1 ) for i in xrange ( m ) : tmp = a [ i ] lapack_routine ( n , n , tmp , n , pivots [ i ] , 0 ) sign = 1. - 2. * ( nm . add . reduce ( pivots != flags , axis = 1 ) % 2 ) idx = nm . arange ( n ) d = a [ : , idx , idx ] absd = nm . absolute ( d ) sign *= nm . multiply . reduce ( d / absd , axis = 1 ) nm . log ( absd , absd ) logdet = nm . add . reduce ( absd , axis = - 1 ) return sign * nm . exp ( logdet ) def get_volume ( el , nd ) : dim = nd . shape [ 1 ] nnd = el . shape [ 1 ] etype = '%d_%d' % ( dim , nnd ) if etype == '2_4' or etype == '3_8' : el = elems_q2t ( el ) nel = el . shape [ 0 ] mul = 1.0 / factorial ( dim ) if dim == 3 : mul *= - 1.0 mtx = nm . ones ( ( nel , dim + 1 , dim + 1 ) , dtype = nm . double ) mtx [ : , : , : - 1 ] = nd [ el , : ] vols = mul * dets_fast ( mtx . copy ( ) ) vol = vols . sum ( ) bc = nm . dot ( vols , mtx . sum ( 1 ) [ : , : - 1 ] / nnd ) bc /= vol return vol , bc n_nod = mesh . n_nod if weights is None : node_group = nm . ones ( ( n_nod , ) , dtype = nm . int8 ) * 2 sndi , edges = get_snodes_uedges ( mesh . conns [ 0 ] , mesh . descs [ 0 ] ) if bconstr : node_group [ sndi ] = 4 end1 = edges [ : , 0 ] end2 = edges [ : , 1 ] idxs = nm . where ( node_group [ end2 ] >= node_group [ end1 ] ) rows1 = end1 [ idxs ] cols1 = end2 [ idxs ] idxs = nm . where ( node_group [ end1 ] >= node_group [ end2 ] ) rows2 = end2 [ idxs ] cols2 = end1 [ idxs ] crows = nm . concatenate ( ( rows1 , rows2 ) ) ccols = nm . concatenate ( ( cols1 , cols2 ) ) costs = sps . coo_matrix ( ( nm . ones_like ( crows ) , ( crows , ccols ) ) , shape = ( n_nod , n_nod ) , dtype = nm . double ) idxs = range ( n_nod ) aux = sps . coo_matrix ( ( 1.0 / nm . asarray ( costs . sum ( 1 ) ) . squeeze ( ) , ( idxs , idxs ) ) , shape = ( n_nod , n_nod ) , dtype = nm . double ) weights = ( aux . tocsc ( ) * costs . tocsc ( ) ) . tocsr ( ) coors = taubin ( mesh . coors , weights , lam , mu , n_iter ) if volume_corr : volume0 , bc = get_volume ( mesh . conns [ 0 ] , mesh . coors ) volume , _ = get_volume ( mesh . conns [ 0 ] , coors ) scale = volume0 / volume coors = ( coors - bc ) * scale + bc return coors
FE mesh smoothing.

    Based on:

    [1] Steven K. Boyd, Ralph Muller, Smooth surface meshing for automated
    finite element model generation from 3D image data, Journal of
    Biomechanics, Volume 39, Issue 7, 2006, Pages 1287-1295,
    ISSN 0021-9290, 10.1016/j.jbiomech.2005.03.006.
    (http://www.sciencedirect.com/science/article/pii/S0021929005001442)

    Parameters
    ----------
    mesh : mesh
        FE mesh.
    n_iter : integer, optional
        Number of iteration steps.
    lam : float, optional
        Smoothing factor, see [1].
    mu : float, optional
        Unshrinking factor, see [1].
    weights : array, optional
        Edge weights, see [1].
    bconstr: logical, optional
        Boundary constraints, if True only surface smoothing performed.
    volume_corr: logical, optional
        Correct volume after smoothing process.

    Returns
    -------
    coors : array
        Coordinates of mesh nodes.