Python pygimli.x函数代码示例

    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