def test_Helmholtz(self):
        """
            d² u / d x² + k u + f = 0
            k = 2
            a) P1(exact)
                u = x
                f = -2x
            b) P2(exact)
                u = x*x
                f = -(2 + 2x*x)
        """
        h = np.pi/2 / 21
        x = np.arange(0.0, np.pi/2, h)
        mesh = pg.createGrid(x)

        ### test a)
        k = 2.0
        u = lambda _x: _x
        f = lambda _x: -(k * u(_x))

        x = pg.x(mesh)
        dirichletBC = [[1, u(min(x))], [2, u(max(x))]]
        uFEM = pg.solve(mesh, a=1, b=k, f=f(x), bc={'Dirichlet': dirichletBC})
        np.testing.assert_allclose(uFEM, u(x))

        ### test b)
        u = lambda _x: _x * _x
        f = lambda _x: -(2. + k *u(_x))

        mesh = mesh.createP2()
        x = pg.x(mesh)
        dirichletBC = [[1, u(min(x))], [2, u(max(x))]]
        uFEM = pg.solve(mesh, a=1, b=k, f=f(x), bc={'Dirichlet': dirichletBC})
        np.testing.assert_allclose(uFEM, u(x), atol=1e-6)

    def ani(i):
        axGra.clear()
        axGra.plot(pg.x(gravPoints), dz[i])
        axGra.plot(pg.x(gravPoints), pg.y(gravPoints), 'v', color='black')
        axGra.set_ylabel('Grav in mGal')
        axGra.set_xlim((-20, 20))
        axGra.set_ylim((0, 0.001))
        axGra.grid()

        pg.mplviewer.setMappableData(gciDDe, abs(dDens[i]),
                                     cMin=0, cMax=20,
                                     logScale=False)

    def showVA(self, ax=None, t=None, name='va', pseudosection=False,
               squeeze=True, full=True):
        """show apparent velocity as image plot

        TODO showXXX commands need to return axes and cbar .. if there is one

        """

        if ax is None:
            fig, ax = plt.subplots()
            self.figs[name] = fig

        self.axs[name] = ax
        if t is None:
            t = self.dataContainer('t')

        px = pg.x(self.dataContainer.sensorPositions())
        gx = np.array([px[int(g)] for g in self.dataContainer("g")])
        sx = np.array([px[int(s)] for s in self.dataContainer("s")])
        offset = self.getOffset(full=full)
        va = offset / t

        if pseudosection:
            midpoint = (gx + sx) / 2
            plotVecMatrix(midpoint, offset, va, squeeze=True, ax=ax,
                          label='Apparent slowness [s/m]')
        else:
            plotVecMatrix(gx, sx, va, squeeze=squeeze, ax=ax,
                          label='Apparent velocity [m/s]')
#        va = showVA(ax, self.dataContainer)
#        plt.show(block=False)
        return va

def cellDataToBoundaryGrad(mesh, v, vGrad):
    """
    """
    if len(v) != mesh.cellCount() or len(vGrad) != mesh.cellCount():
        raise
    gB = mesh.cellDataToBoundaryGradient(v, vGrad)
    return np.vstack([pg.x(gB), pg.y(gB), pg.z(gB)]).T

    gB = np.zeros((mesh.boundaryCount(), 3))

    for b in mesh.boundaries():
        leftCell = b.leftCell()
        rightCell = b.rightCell()
        gr = pg.RVector3(0.0, 0.0, 0.0)
        t = (b.node(1).pos() - b.node(0).pos()).norm()

        if leftCell and rightCell:
            df1 = b.center().distance(leftCell.center())
            df2 = b.center().distance(rightCell.center())

            gr = b.norm() * \
                (v[rightCell.id()] - v[leftCell.id()]) / (df1 + df2)

            grL = t * t.dot(vGrad[leftCell.id()])
            grR = t * t.dot(vGrad[rightCell.id()])

            gr += (grL + grR) * 0.5

        elif leftCell:
            gr = t * t.dot(vGrad[leftCell.id()])

        gB[b.id(), 0] = gr[0]
        gB[b.id(), 1] = gr[1]
        gB[b.id(), 2] = gr[2]
    return gB

def cellDataToBoundaryGrad(mesh, v, vGrad):
    """TODO Documentme."""
    if len(v) != mesh.cellCount() or len(vGrad) != mesh.cellCount():
        raise BaseException("len(v) dismatch mesh.cellCount()")

    gB = mesh.cellDataToBoundaryGradient(v, vGrad)
    return np.vstack([pg.x(gB), pg.y(gB), pg.z(gB)]).T

def cellDataToCellGrad(mesh, v, CtB):
    """TODO Documentme."""
    if len(v) != mesh.cellCount():
        print(len(v), mesh.cellCount())
        raise BaseException("len of v missmatch mesh.cellCount()")
    div = mesh.boundaryDataToCellGradient(CtB * v)
    return np.vstack([pg.x(div), pg.y(div), pg.z(div)]).T

def createCoarsePoly( coarseData ):
    boundary = 1250.0
    mesh = g.Mesh()

    x = g.x( coarseData )
    y = g.y( coarseData )
    z = g.z( coarseData )

    xMin = min( x ); xMax = max( x )    
    yMin = min( y ); yMax = max( y )
    zMin = min( z ); zMax = max( z )        

    print(xMin, xMax, yMin, yMax)
    border = max( (xMax - xMin) * boundary / 100.0, (yMax - yMin) * boundary / 100.0);

    n1 = mesh.createNode( xMin - border, yMin - border, zMin, 1 )
    n2 = mesh.createNode( xMax + border, yMin - border, zMin, 2 )
    n3 = mesh.createNode( xMax + border, yMax + border, zMin, 3 )
    n4 = mesh.createNode( xMin - border, yMax + border, zMin, 4 )
  
    mesh.createEdge( n1, n2, 12 );
    mesh.createEdge( n2, n3, 23 );
    mesh.createEdge( n3, n4, 34 );
    mesh.createEdge( n4, n1, 41 );

    for p in coarseData:
        mesh.createNode( p )

    return mesh

def plotFirstPicks(ax, data, tt=None, plotva=False, marker='x-'):
    """plot first arrivals as lines"""
    px = pg.x(data.sensorPositions())
    gx = np.array([px[int(g)] for g in data("g")])
    sx = np.array([px[int(s)] for s in data("s")])
    if tt is None:
        tt = np.array(data("t"))
    if plotva:
        tt = np.absolute(gx - sx) / tt

    uns = np.unique(sx)

    cols = plt.cm.tab10(np.arange(10))

    for i, si in enumerate(uns):
        ti = tt[sx == si]
        gi = gx[sx == si]
        ii = gi.argsort()
        ax.plot(gi[ii], ti[ii], marker, color=cols[i % 10])
        ax.plot(si, 0., 's', color=cols[i % 10], markersize=8)

    ax.grid(True)
    if plotva:
        ax.set_ylabel("Apparent velocity (m/s)")
    else:
        ax.set_ylabel("Traveltime (s)")
    ax.set_xlabel("x (m)")
    ax.invert_yaxis()

def integrate(f, ent, order):
    """ integrate function """

    J = 0
    x = []
    w = []
    if type(ent) is list:

        a = ent[0]
        b = ent[1]
        xs = pg.IntegrationRules.instance().gauAbscissa(order)
        w = pg.IntegrationRules.instance().gauWeights(order)

        x = (b - a) / 2.0 * pg.x(xs) + (a + b) / 2.0
        J = (b - a) / 2.0

    else:
        J = ent.shape().jacobianDeterminant()
        xs = pg.IntegrationRules.instance().abscissa(ent.shape(), order)
        w = pg.IntegrationRules.instance().weights(ent.shape(), order)

        x = [ent.shape().xyz(xsi) for xsi in xs]

    for xx in x:
        print xx
    print w
    print funct(f)(x)
    return J * sum(funct(f)(x) * w)

def create_mesh_and_data(n):
    nc = np.linspace(-2.0, 0.0, n)
    mesh = pg.createMesh2D(nc, nc)
    mcx = pg.x(mesh.cellCenter())
    mcy = pg.y(mesh.cellCenter())
    data = np.cos(1.5 * mcx) * np.sin(1.5 * mcy)
    return mesh, data

def cellDataToCellGrad(mesh, v, CtB):
    if len(v) != mesh.cellCount():
        print(len(v), mesh.cellCount())
        raise
    div = mesh.boundaryDataToCellGradient(CtB * v)
    return np.vstack([pg.x(div), pg.y(div), pg.z(div)]).T

    vF = cellDataToBoundaryData(mesh, v)
    gC = np.zeros((mesh.cellCount(), 3))

    for b in mesh.boundaries():

        leftCell = b.leftCell()
        rightCell = b.rightCell()
        vec = b.norm() * vF[b.id()] * b.size()

        if leftCell:
            gC[leftCell.id(), 0] += vec[0]
            gC[leftCell.id(), 1] += vec[1]
            gC[leftCell.id(), 2] += vec[2]
        if rightCell:
            gC[rightCell.id(), 0] -= vec[0]
            gC[rightCell.id(), 1] -= vec[1]
            gC[rightCell.id(), 2] -= vec[2]

    gC[:, 0] /= mesh.cellSizes()
    gC[:, 1] /= mesh.cellSizes()
    gC[:, 2] /= mesh.cellSizes()

    return gC

def transform2DMeshTo3D(mesh, x, y, z=None):
    """
    Transform a 2D mesh into 3D coordinates using a point list (e.g. from GPS)

    Parameters
    ----------
    mesh: GIMLi::Mesh
    x,y: array of x/y positions along 2d profile
    z: optional height to add (topographical correction if computed flat earth)

    See Also
    --------

    References
    ----------
    """

    # get mesh node positions
    mt, mz = pg.x( mesh.positions() ), pg.y( mesh.positions() ) # mesh tape and z
    # compute length of reference points along tape
    pt = np.hstack( (0., np.cumsum( np.sqrt( np.diff( x )**2 + np.diff( y )**2 ) ) ) )
    #  interpolate node positions from tape to x/y using tape positions
    mx = np.interp( mt, pt, x )
    my = np.interp( mt, pt, y )
    # compute z offset by interpolating z
    if z is None:
        oz = np.zeros( len(mt) )
    else:
        oz = np.interp( mt, pt, z )

    # set the positions in the mesh
    for i, node in enumerate( mesh.nodes() ):
        node.setPos( pg.RVector3( mx[i], my[i], mz[i]+oz[i] ) )

    def test_createPartMesh(self):
        mesh = pg.createMesh1D(np.linspace(0, 1, 10))
        self.assertEqual(mesh.cellCount(), 9)

        mesh2 = mesh.createMeshByCellIdx(
            pg.find(pg.x(mesh.cellCenters()) < 0.5))
        self.assertEqual(mesh2.cellCount(), 4)
        self.assertEqual(mesh2.cellCenters()[-1][0] < 0.5, True)

    def getMidpoint(self, data=None):
        """Return vector of offsets (in m) between shot and receiver."""
        if data is None:
            data = self.dataContainer

        px = pg.x(data.sensorPositions())
        gx = np.array([px[int(g)] for g in data("g")])
        sx = np.array([px[int(s)] for s in data("s")])
        return (gx + sx) / 2

        def _testP1_(mesh, show=False):
            """ Laplace u = 0 solves u = x for u(r=0)=0 and u(r=1)=1
                Test for u == exact x for P1 base functions
            """
            u = pg.solve(mesh, a=1, b=0, f=0, 
                         bc={'Dirichlet': [[1, 0], [2, 1]]})

            if show:
                if mesh.dim()==1:    
                    pg.plt.plot(pg.x(mesh), u)
                    pg.wait()
                elif mesh.dim()==2:
                    pg.show(mesh, u, label='u')
                    pg.wait()

            xMin = mesh.xmin()
            xSpan = (mesh.xmax() - xMin)
            np.testing.assert_allclose(u, (pg.x(mesh)-xMin)/ xSpan)
            return u

def calcInvBlock(mesh, dens, out='gravInv'):

    # extract block delta density
    densBlock = pg.RVector(dens)
    densMarker2 = dens[pg.find(mesh.cellMarker() == 2)[0]]
#    densBlock[(mesh.cellMarker() == 1)|(mesh.cellMarker() == 3)] = densMarker2
    densBlock[pg.find((mesh.cellMarker() == 1) | (mesh.cellMarker() == 3))] = \
        densMarker2
    densBlock -= densMarker2

    # define meausrement positions
    gravPointsX = np.linspace(-20, 20, 41)
    sensorPositions = np.vstack((gravPointsX, np.zeros(len(gravPointsX)))).T

    # solve analytical
    gz = solveGravimetry(mesh, densBlock, pnts=sensorPositions, complete=False)

    # noisyfy
    errAbs = 0.00001
    dzerr = np.random.randn(len(sensorPositions)) * errAbs
    gz = gz + dzerr

    # createParamesh
    paraMesh = pg.createGrid(x=np.linspace(-20, 20, 41),
                             y=np.linspace(-20, 0, 21))

    # init Gravimetry manager (should do meshing, simulation and noisying)
    Grav = Gravimetry(verbose=True)

    model = Grav.invert(sensorPositions, gz, errAbs, verbose=1, mesh=paraMesh)

    fig, ax = plt.subplots()
    ax.plot(pg.x(sensorPositions), gz, label='gz')
    ax.plot(pg.x(sensorPositions), Grav.inv.response(), label='response')
    ax.legend()
    ax.grid()
    ax.set_xlabel('$x$ [m]')
    ax.set_ylabel('$\partial u / \partial z$ [mGal]')
    plt.show(block=False)
    ax.figure.savefig(out, bbox_inches='tight')

    return Grav, densBlock

def invertGravimetry(gravPoints, dz):

    dzerr = np.random.randn(len(gravPoints)) * 0.0001
    dz = dz + dzerr

    mesh = pg.createGrid(x=np.linspace(-20, 20, 41),
                         y=np.linspace(-20, 0, 21))

    grav = Gravimetry(verbose=True)

    model = grav.invert(gravPoints, dz, verbose=1, mesh=mesh)

    plt.plot(pg.x(gravPoints), dz)
    plt.plot(pg.x(gravPoints), grav.inv.response())

    paraDomain=grav.fop.regionManager().paraDomain()
    pg.show(paraDomain, model, colorBar=1, hold=1)

    pg.showNow()
    plt.show()
    pass

def createGradientModel2D(data, mesh, VTop, VBot):
    """
    Create 2D velocity gradient model.

    Creates a smooth, linear, starting model that takes the slope
    of the topography into account. This is done by fitting a straight line
    and using the distance to that as the depth value.
    Known as "The Marcus method"


    Parameters
    ----------
    data : pygimli DataContainer
        The topography list is in here.
    mesh : pygimli.Mesh
        The parametric mesh used for the inversion
    VTop : float
        The velocity at the surface of the mesh
    VBot : float
        The velocity at the bottom of the mesh

    Returns
    -------
    model : pygimli Vector, length M
        A numpy array with slowness values that can be used to start
        the inversion.
    """

    p = np.polyfit(pg.x(data.sensorPositions()), pg.y(data.sensorPositions()),
                   deg=1)  # slope-intercept form
    n = np.asarray([-p[0], 1.0])  # normal vector
    nLen = np.sqrt(np.dot(n, n))

    x = pg.x(mesh.cellCenters())
    z = pg.y(mesh.cellCenters())
    pos = np.column_stack((x, z))
    d = np.array([np.abs(np.dot(pos[i, :], n) - p[1]) / nLen
                  for i in range(pos.shape[0])])

    return np.interp(d, [min(d), max(d)], [1.0 / VTop, 1.0 / VBot])

 def getOffset(self, full=False):
     """return vector of offsets (in m) between shot and receiver"""
     if full:
         pos = self.dataContainer.sensorPositions()
         s, g = self.dataContainer('s'), self.dataContainer('g')
         nd = self.dataContainer.size()
         off = [pos[int(s[i])].distance(pos[int(g[i])]) for i in range(nd)]
         return np.absolute(off)
     else:
         px = pg.x(self.dataContainer.sensorPositions())
         gx = np.array([px[int(g)] for g in self.dataContainer("g")])
         sx = np.array([px[int(s)] for s in self.dataContainer("s")])
         return np.absolute(gx - sx)

def drawTravelTimeData(axes, data, t=None):
    """
        Draw first arrival traveltime data into mpl axes a.
        data of type \ref DataContainer must contain sensorIdx 's' and 'g'
        and thus being numbered internally [0..n)
    """

    x = pg.x(data.sensorPositions())
#    z = pg.z(data.sensorPositions())

    shots = pg.unique(pg.sort(data('s')))
    geoph = pg.unique(pg.sort(data('g')))

    startOffsetIDX = 0

    if min(min(shots), min(geoph)) == 1:
        startOffsetIDX = 1

    tShow = data('t')
    if t is not None:
        tShow = t

    axes.set_xlim([min(x), max(x)])
    axes.set_ylim([max(tShow), -0.002])
    axes.figure.show()

    for shot in shots:
        gIdx = pg.find(data('s') == shot)
        sensorIdx = [int(i__ - startOffsetIDX) for i__ in data('g')[gIdx]]
        axes.plot(x[sensorIdx], tShow[gIdx], 'x-')

    yPixel = axes.transData.inverted().transform_point((1, 1))[1] - \
        axes.transData.inverted().transform_point((0, 0))[1]
    xPixel = axes.transData.inverted().transform_point((1, 1))[0] - \
        axes.transData.inverted().transform_point((0, 0))[0]

    # draw shot points
    axes.plot(x[[int(i__ - startOffsetIDX) for i__ in shots]],
           np.zeros(len(shots)) + 8. * yPixel, 'gv', markersize=8)

    # draw geophone points
    axes.plot(x[[int(i__ - startOffsetIDX) for i__ in geoph]],
           np.zeros(len(geoph)) + 3. * yPixel, 'r^', markersize=8)

    axes.grid()
    axes.set_ylim([max(tShow), +16. * yPixel])
    axes.set_xlim([min(x) - 5. * xPixel, max(x) + 5. * xPixel])

    axes.set_xlabel('x-Coordinate [m]')
    axes.set_ylabel('Traveltime [ms]')