基于Python的二维谱元法声波方程正演计算

谱元法其实一种高阶的有限元方法，本质上还是属于伽辽金方法的一种，首先按照有限单元的思想对空间进行离散，然后再单元内部进行插值。同时，利用某些配点和积分方案完成区域的积分。

谱元法与有限元法的区别和联系？

首先有限元的单元形式多种多样，而谱元法一般适用在二维四边形和三维六面体中，不能像有限元方法的三角单元那样。另外，谱元法选取的插值基函数为拉格朗日基函数，有限元采用的是线性插值基函数。在单元内部，谱元法使用gauss lobatto legendre（简称gll）积分方案，将质量矩阵和刚度矩阵变为一系列权重系数的组合，同时由于多项式的正交性，最终的质量矩阵为对角矩阵。

1.1 谱元法求解声波方程的步骤

1. 声波方程的弱形式

$$\begin{array}{rl}{\partial }_t^{2}u& =v^2{\nabla }^{2}u+f\\ \iint {\partial }_t^{2}u\cdot \alpha d\Omega & =\iint v^2{\nabla }^{2}u\cdot \alpha d\Omega +\iint f\cdot \alpha d\Omega \end{array}$$其中$\alpha (x,y)$为测试函数

2. 使用分部积分，并应用自由边界条件

$$\begin{array}{rl}\iint v^2{\nabla }^{2}u\cdot \alpha d\Omega & ={\left[v^2\nabla u\cdot \alpha \right]|}_{l_1}^{l_2}-v^2\iint \nabla u\nabla \alpha d\Omega \\ & =-v^2\iint \nabla u\nabla \alpha d\Omega ({\left[v^2\nabla u\cdot \alpha \right]|}_{l_1}^{l_2}=0)\end{array}$$

3. 对解进行插值近似

$$\begin{array}{rl}& u(t,x,y)=\sum _{i=0}^N{u}_i(t){\ell }_i(x,y)\\ \iint {\partial }_t^2\left[\sum _{i=0}^N{u}_i(t){\ell }_i(x,y)\right]\cdot \alpha d\Omega & +v^2\iint \nabla \left[\sum _{i=0}^N{u}_i(t){\ell }_i(x,y)\right]\cdot \nabla \alpha d\Omega =\iint f(t,x,y)\cdot \alpha d\Omega \\ \sum _{i=0}^N{\partial }_t^2{u}_i(t)\iint {\ell }_i(x,y)\alpha (x,y)d\Omega & +v^2\sum _{i=0}^N{u}_i(t)\iint \nabla {\ell }_i(x,y)\nabla \alpha (x,y)d\Omega =\iint f(t,x,y)\cdot \alpha d\Omega \end{array}$$其中${\ell }_i(x,y)={\ell }_i(x)\times {\ell }_i(y)$为二维的拉格朗日插值基函数。

4. 应用伽辽金法，令测试函数为相同的基函数

$$\begin{array}{rl}\sum _{i=0}^N{\partial }_t^2{u}_i(t)\iint {\ell }_i(x,y){\ell }_j(x,y)d\Omega & +v^2\sum _{i=0}^N{u}_i(t)\iint \nabla {\ell }_i(x,y)\nabla {\ell }_j(x,y)d\Omega =\iint f(t,x,y){\ell }_j(x,y)d\Omega \end{array}$$

5. 将原方程写成弱形式的矩阵方程组

$$\begin{array}{rl}& {\mathbf{M}}^T{\partial }_t^2\mathbf{u}+{\mathbf{K}}^T\mathbf{u}=\mathbf{f}\\ {\mathbf{M}}^T& =\iint {\ell }_i(x,y){\ell }_j(x,y)d\Omega \\ {\mathbf{K}}^T& =\iint \nabla {\ell }_i(x,y)\nabla {\ell }_j(x,y)d\Omega \\ \mathbf{f}& =\iint f(t,x,y){\ell }_j(x,y)d\Omega \end{array}$$

6. 使用二阶中心差分对时间进行延拓

$$\begin{array}{r}\mathbf{u}(t+dt)=d{t}^2({\mathbf{M}}^T{)}^{-1}\left[\mathbf{f}-{\mathbf{K}}^T\mathbf{u}\right]+2\mathbf{u}-\mathbf{u}(t-dt)\end{array}$$

1.2 拉格朗日多项式与数值积分

根据数值定理，对于 $N+1$ 个点有不超过 $N$ 次多项式与之对应，拉格朗日多项式满足该条件，且具有以下形式：
$$\begin{array}{r}{\ell }_i^N(x)=\prod _{i\ne j}^{N+1}\frac{x-x_j}{x_i-x_j},i,j=1,2,\dots ,N+1\end{array}$$其特点是在当前节点处值为1，其余节点处为零。
连续函数的积分可以用可以解析积分的多项式近似代替来计算。我们再次使用拉格朗日多项式作为插值函数，得到高斯-洛巴托-勒让德积分。在这里，gll点被用来计算积分。
$$\begin{array}{r}{\int }_{-1}^{1}f(x)dx\approx {\int }_{-1}^1{P}_N(x)dx=\sum _{i=1}^{N+1}{\omega }_{i}f(x_i)\end{array}$$其中$x_i$为gll点，${\omega }_i$为相应gll点的权重系数。

1.3 网格划分

本文使用的网格是长宽相同的正方形单元网格，单元内部使用gll点进行划分。

网格划分

正方形单元映射

1.4 质量矩阵的计算

写出质量矩阵的计算式：
$$\begin{array}{rl}M& =\iint {\ell }_i(x,y){\ell }_j(x,y)d\Omega \\ & ={\int }_{-1}^1{\int }_{-1}^1{\ell }_i(u,v){\ell }_j(u,v)\left|J\right|dudv\\ & ={\int }_{-1}^1\sum _{k=0}^N{\omega }_i{\ell }_i(u_k,v){\ell }_j(u_k,v)\left|J\right|dv\\ & =\sum _{k=0}^N\sum _{l=0}^N{\omega }_i{\omega }_j{\ell }_i(u_k,v_l){\ell }_j(u_k,v_l)\left|J\right|\\ & =\sum _{k=0}^N\sum _{l=0}^N{\omega }_i{\omega }_j{\ell }_i(u_k){\ell }_i(v_l){\ell }_j(u_k){\ell }_j(v_l)\left|J\right|\\ & ={\omega }_i{\omega }_j{\delta }_{ij}\left|J\right|\end{array}$$由此可知，谱元法的质量矩阵是一个对角矩阵，且元素为对应位置的积分权重系数的乘积。由于本文使用的网格中每一个单元形状一样，所以都有相同的局部质量矩阵，计算可得该局部单元的质量矩阵为：

局部质量矩阵及对角元素

1.5 刚度矩阵的计算

写出刚度矩阵的计算式：
$$\begin{array}{rl}K& =\iint \nabla {\ell }_i(x,y)\nabla {\ell }_j(x,y)d\Omega \\ & ={\int }_{-1}^1{\int }_{-1}^1\nabla {\ell }_i(x,y)\nabla {\ell }_j(x,y)\left|J\right|dudv\\ & ={\int }_{-1}^1{\int }_{-1}^1\left[\nabla {\ell }_i(u,v)J^{-1}\right]\left[\nabla {\ell }_j(u,v)J^{-1}\right]\left|J\right|dudv\\ & ={\int }_{-1}^1{\int }_{-1}^1\left\{\nabla \left[{\ell }_i(u){\ell }_i(v)\right]J^{-1}\right\}\left\{\nabla \left[{\ell }_j(u){\ell }_j(v)\right]J^{-1}\right\}\left|J\right|dudv\\ & ={\int }_{-1}^1{\int }_{-1}^1\left\{\left[{\partial }_u{\ell }_i(u){\ell }_i(v),{\partial }_v{\ell }_i(u){\ell }_i(v)\right]J^{-1}\right\}\left\{\left[{\partial }_u{\ell }_j(u){\ell }_j(v),{\partial }_v{\ell }_j(u){\ell }_j(v)\right]J^{-1}\right\}\left|J\right|dudv\\ & ={\int }_{-1}^1{\int }_{-1}^1\frac{2}{le}\left[{\partial }_u{\ell }_i(u){\ell }_i(v),{\partial }_v{\ell }_i(u){\ell }_i(v)\right]\frac{2}{le}{\left[{\partial }_u{\ell }_j(u){\ell }_j(v),{\partial }_v{\ell }_j(u){\ell }_j(v)\right]}^T\left|J\right|dudv(J^{-1}=\left[\begin{array}{cc}\frac{2}{le}& 0\\ 0& \frac{2}{le}\end{array}\right])\\ & ={\int }_{-1}^1{\int }_{-1}^1\left[{\partial }_u{\ell }_i(u){\ell }_i(v),{\partial }_v{\ell }_i(u){\ell }_i(v)\right]{\left[{\partial }_u{\ell }_j(u){\ell }_j(v),{\partial }_v{\ell }_j(u){\ell }_j(v)\right]}^{T}dudv(\left|J\right|=(\frac{le}{2}{)}^2)\\ & ={\int }_{-1}^1\sum _{k=0}^N{\omega }_k\left[{\partial }_u{\ell }_i(u_k){\ell }_i(v),{\partial }_v{\ell }_i(u_k){\ell }_i(v)\right]{\left[{\partial }_u{\ell }_j(u_k){\ell }_j(v),{\partial }_v{\ell }_j(u_k){\ell }_j(v)\right]}^{T}dv\\ & =\sum _{k=0}^N\sum _{l=0}^N{\omega }_k{\omega }_l\left[{\partial }_u{\ell }_i(u_k){\ell }_i(v_l),{\partial }_v{\ell }_i(u_k){\ell }_i(v_l)\right]{\left[{\partial }_u{\ell }_j(u_k){\ell }_j(v_l),{\partial }_v{\ell }_j(u_k){\ell }_j(v_l)\right]}^T\\ & =\sum _{k=0}^N\sum _{l=0}^N{\omega }_k{\omega }_l\left[{\partial }_u{\ell }_i(u_k){\ell }_i(v_l){\partial }_u{\ell }_j(u_k){\ell }_j(v_l)+{\partial }_v{\ell }_i(u_k){\ell }_i(v_l){\partial }_v{\ell }_j(u_k){\ell }_j(v_l)\right]\\ & =\sum _{k=0}^N\sum _{l=0}^N{\omega }_k{\omega }_l\left[{\partial }_u{\ell }_i(u_k){\partial }_u{\ell }_j(u_k){\delta }_{iu,ju}+{\partial }_v{\ell }_i(u_k){\partial }_v{\ell }_j(u_k){\delta }_{iv,jv}\right]\end{array}$$
即局部刚度矩阵计算只在$i,j$ 点横纵坐标至少有一个相等的情况下才有值，一个都不相等则为零。所以局部刚度计算可按三种情况进行计算：

局部刚度矩阵非零的三种情况 由此，我们可以计算得到局部刚度矩阵：

二、数值实验

2.1 参数设置

本次实验计算的为均匀模型下的弹性波（与声波无异），具体参数设置如下：

2.2 总装质量矩阵与刚度矩阵

将之前计算的局部质量矩阵与局部刚度矩阵按照单元进行组装，得到全局的质量矩阵与刚度矩阵如下（绘制部分）：

全局矩阵

2.3 震源

本次实验采用高斯震源，具有以下形式：

高斯震源

2.4 实验结果

将计算好的矩阵进行显示格式迭代，可得各个时刻的波场值：

高斯震源

三、示例代码

#!/usr/bin/env python
# -*- coding: utf-8 -*-
'''
作者：@fm
'''


def gll(N):
    """
    Returns GLL (Gauss Lobato Legendre module with collocation points and
    weights)
    """
    # Initialization of integration weights and collocation points
    # [xi, weights] =  gll(N)
    # Values taken from Diploma Thesis Bernhard Schuberth
    if N == 2:
        xi = [-1.0, 0.0, 1.0]
        weights = [0.33333333, 1.33333333, 0.33333333]
    elif N == 3:
        xi = [-1.0, -0.447213595499957, 0.447213595499957, 1.0]
        weights = [0.1666666667, 0.833333333, 0.833333333, 0.1666666666]
    elif N == 4:
        xi = [-1.0, -0.6546536707079772, 0.0, 0.6546536707079772, 1.0]
        weights = [0.1, 0.544444444, 0.711111111, 0.544444444, 0.1]
    elif N == 5:
        xi = [-1.0, -0.7650553239294647, -0.285231516480645, 0.285231516480645,
              0.7650553239294647, 1.0]
        weights = [0.0666666666666667,  0.3784749562978470,
                   0.5548583770354862, 0.5548583770354862, 0.3784749562978470,
                   0.0666666666666667]
    elif N == 6:
        xi = [-1.0, -0.8302238962785670, -0.4688487934707142, 0.0,
              0.4688487934707142, 0.8302238962785670, 1.0]
        weights = [0.0476190476190476, 0.2768260473615659, 0.4317453812098627,
                   0.4876190476190476, 0.4317453812098627, 0.2768260473615659,
                   0.0476190476190476]
    elif N == 7:
        xi = [-1.0, -0.8717401485096066, -0.5917001814331423,
              -0.2092992179024789, 0.2092992179024789, 0.5917001814331423,
              0.8717401485096066, 1.0]
        weights = [0.0357142857142857, 0.2107042271435061, 0.3411226924835044,
                   0.4124587946587038, 0.4124587946587038, 0.3411226924835044,
                   0.2107042271435061, 0.0357142857142857]
    elif N == 8:
        xi = [-1.0, -0.8997579954114602, -0.6771862795107377,
              -0.3631174638261782, 0.0, 0.3631174638261782,
              0.6771862795107377, 0.8997579954114602, 1.0]
        weights = [0.0277777777777778, 0.1654953615608055, 0.2745387125001617,
                   0.3464285109730463, 0.3715192743764172, 0.3464285109730463,
                   0.2745387125001617, 0.1654953615608055, 0.0277777777777778]
    elif N == 9:
        xi = [-1.0, -0.9195339081664589, -0.7387738651055050,
              -0.4779249498104445, -0.1652789576663870, 0.1652789576663870,
              0.4779249498104445, 0.7387738651055050, 0.9195339081664589, 1.0]
        weights = [0.0222222222222222, 0.1333059908510701, 0.2248893420631264,
                   0.2920426836796838, 0.3275397611838976, 0.3275397611838976,
                   0.2920426836796838, 0.2248893420631264, 0.1333059908510701,
                   0.0222222222222222]
    elif N == 10:
        xi = [-1.0, -0.9340014304080592, -0.7844834736631444,
              -0.5652353269962050, -0.2957581355869394, 0.0,
              0.2957581355869394, 0.5652353269962050, 0.7844834736631444,
              0.9340014304080592, 1.0]
        weights = [0.0181818181818182, 0.1096122732669949, 0.1871698817803052,
                   0.2480481042640284, 0.2868791247790080, 0.3002175954556907,
                   0.2868791247790080, 0.2480481042640284, 0.1871698817803052,
                   0.1096122732669949, 0.0181818181818182]
    elif N == 11:
        xi = [-1.0, -0.9448992722228822, -0.8192793216440067,
              -0.6328761530318606, -0.3995309409653489, -0.1365529328549276,
              0.1365529328549276, 0.3995309409653489, 0.6328761530318606,
              0.8192793216440067, 0.9448992722228822, 1.0]
        weights = [0.0151515151515152, 0.0916845174131962, 0.1579747055643701,
                   0.2125084177610211, 0.2512756031992013, 0.2714052409106962,
                   0.2714052409106962, 0.2512756031992013, 0.2125084177610211,
                   0.1579747055643701, 0.0916845174131962, 0.0151515151515152]
    elif N == 12:
        xi = [-1.0, -0.9533098466421639, -0.8463475646518723,
              -0.6861884690817575, -0.4829098210913362, -0.2492869301062400,
              0.0, 0.2492869301062400, 0.4829098210913362,
              0.6861884690817575, 0.8463475646518723, 0.9533098466421639,
              1.0]
        weights = [0.0128205128205128, 0.0778016867468189, 0.1349819266896083,
                   0.1836468652035501, 0.2207677935661101, 0.2440157903066763,
                   0.2519308493334467, 0.2440157903066763, 0.2207677935661101,
                   0.1836468652035501, 0.1349819266896083, 0.0778016867468189,
                   0.0128205128205128]
    else:
        raise NotImplementedError

    return xi, weights

# -*- coding: utf-8 -*-
'''
作者：@fm
'''
import numpy as np
from gll import gll

def creat_mesh(le, nx, ny, N):
    '''
    element的大小
    nx:x方向上的element的数量
    ny:y方向上的element的数量
    N:拉格朗日多项式的阶数
    用于谱元法的网格剖分
    '''
    [xi,w]=gll(N=N)
    for i in range(len(xi)):
        xi[i]=(xi[i]+1)/2*le
    # node的数量
    numN = (N*nx+1)*(N*ny+1)
    # element的数量
    numE = nx*ny
    # 一个element有多少nodes
    NofE = (N+1)*(N+1)
    # 二维坐标
    D = 2
    # nodes 坐标
    NC = np.zeros([numN, D])
    
    # nodes 坐标计算
    n = 0
    for i in range(0, N*ny+1):
        for j in range(0, N*nx+1):
            NC[n, 0] = xi[j%N]+(j//N)*le
            NC[n, 1] = xi[i%N]+(i//N)*le

            n += 1
    # element 索引，一个element所有的节点进行索引
    n2 = 0
    EI = np.zeros([numE,(N+1),(N+1)])
    for i in range(ny):
        for j in range(nx):
            for ix in range(N+1):
                
                EI[n2,ix,0:(N+1)]=np.reshape(np.array(range(1,N+2)),[N+1,])+j*N+i*N*(N*nx+1)+ix*(N*nx+1)
            n2+=1
                
    EI = EI.astype(int)    
    return NC, EI

#!/usr/bin/env python
# -*- coding: utf-8 -*-
'''
作者：@fm
'''
import numpy as np

from gll import gll
from lagrange import lagrange
from legendre import legendre


def lagrange1st(N):
    """
    # Calculation of 1st derivatives of Lagrange polynomials
    # at GLL collocation points
    # out = legendre1st(N)
    # out is a matrix with columns -> GLL nodes
    #                        rows -> order
    """
    out = np.zeros([N+1, N+1])

    [xi, w] = gll(N)

    # initialize dij matrix (see Funaro 1993)
    d = np.zeros([N + 1, N + 1])

    for i in range(-1, N):
        for j in range(-1, N):
            if i != j:
                d[i + 1, j + 1] = legendre(N, xi[i + 1]) / \
                    legendre(N, xi[j + 1]) * 1.0 / (xi[i + 1] - xi[j + 1])
            if i == -1:
                if j == -1:
                    d[i + 1, j + 1] = -1.0 / 4.0 * N * (N + 1)
            if i == N-1:
                if j == N-1:
                    d[i + 1, j + 1] = 1.0 / 4.0 * N * (N + 1)

    # Calculate matrix with 1st derivatives of Lagrange polynomials
    for n in range(-1, N):
        for i in range(-1, N):
            sum = 0
            for j in range(-1, N):
                sum = sum + d[i + 1, j + 1] * lagrange(N, n, xi[j + 1])

            out[n + 1, i + 1] = sum
    return(out)

#!/usr/bin/env python
# -*- coding: utf-8 -*-
'''
作者：@fm
'''


def lagrange(N, i, x):
    """
    Function to calculate  Lagrange polynomial for order N and polynomial
    i[0, N] at location x.
    """
    from gll import gll

    [xi, weights] = gll(N)
    fac = 1
    for j in range(-1, N):
        if j != i:
            fac = fac * ((x - xi[j + 1]) / (xi[i + 1] - xi[j + 1]))
    return fac

# -*- coding: utf-8 -*-
'''
作者：@fm
'''
import numpy as np
import matplotlib.pyplot as plt 
import os
import shutil
from gll import gll
from lagrange1st import lagrange1st 
from ricker import ricker
from creat_grid_sem import creat_mesh
from lagrange import lagrange
filepath='images_i21'
nt    = 1000         # number of time steps
xmax  = 600.        # Length of domain [m]
vs    = 2500.         # S velocity [m/s]
rho   = 2000          # Density [kg/m^3]
mu    = rho * vs**2   # Shear modulus mu
N     = 4           # Order of Lagrange polynomials
nx    = 30       # Number of x elements
ny    = 30           # Number of y elements
Tdom  = .1            # Dominant period of Ricker source wavelet
iplot = 30            # Plotting each iplot snapshot
# Space domain
le = xmax/nx        # Length of elements
NC, EI =  creat_mesh(le=le,nx=nx,ny=ny,N=N)
numN = np.size(NC,0)
# 一个单元有多少个节点
nodes_of_one_ele = np.size(EI,1)*np.size(EI,2)
# 计算最小空间间距
dx_min = min(np.diff(NC[0:N+1,0]))
dx_max = max(np.diff(NC[0:N+1,0]))
print(dx_min,dx_max)
# stability limit
eps = 0.1
dt = eps*dx_min/vs
# jacobian map
DJ = (le/2)**2 # (le/2)**2
[xi, wi] = gll(N)
[xj, wj] = gll(N)  
print('dt={:.4f},nt={:.4f}'.format(dt,nt))
def RemoveDir(filepath):
    '''
    如果文件夹不存在就创建，如果文件存在就清空！
    
    '''
    if not os.path.exists(filepath):
        os.mkdir(filepath)
    else:
        shutil.rmtree(filepath)
        os.mkdir(filepath)
def creat_mass_matrix(NC,EI):
    '''
    计算质量矩阵的逆
    '''
    EI =EI -1
    # 节点shul
    
    # 局部质量矩阵,只存储对角元素
    Me = np.zeros(nodes_of_one_ele, dtype =  float)
    n = 0
    for j in range(np.size(EI,2)):
        for i in range(np.size(EI,1)):
            Me[n]=rho*wi[i]*wj[j]*DJ
            n +=1
    M = np.zeros(numN, dtype =  float)
    # 组装质量矩阵
    for ele in EI:
        M[ele.flatten()] += Me
    Minv = np.identity(numN, dtype =  float)
    # 计算质量矩阵的逆
    for i in range(numN):
        Minv[i,i]=1./M[i]
    return Minv
        
        
            
    
def creat_stiffness_matrix(NC,EI):
    '''
    建立刚度矩阵
    '''
    EI=EI-1
    # 节点数量
    numN = np.size(NC,0)
    # 一个单元具有多少节点
    nodes_of_one_ele=np.size(EI,1)*np.size(EI,2)
    l1df = lagrange1st(N)
    # 全局刚度矩阵
    K = np.zeros([numN,numN], dtype =  float)
    # 局部刚度矩阵
    Ke = np.zeros([nodes_of_one_ele,nodes_of_one_ele], dtype =  float)
    # 计算局部刚度矩阵
    for i in range (nodes_of_one_ele):
        for j in range(nodes_of_one_ele):
            ix = i%(N+1)
            iy = i//(N+1)
            
            jx = j%(N+1)
            jy = j//(N+1)
            # 计算三种非零情况
            if ix==jx and iy==jy:
                for k in range(np.size(EI,1)):
                    Ke[i,j]+= mu*wi[k]*wj[jy]*l1df[ix,k]*l1df[jx,k]
                for l in range(np.size(EI,2)):
                    Ke[i,j]+= mu*wi[ix]*wj[l]*l1df[iy,l]*l1df[jy,l]
            elif ix==jx and not iy==jy:
                
                for l in range(np.size(EI,2)):
                    Ke[i,j]+= mu*wi[ix]*wj[l]*l1df[iy,l]*l1df[jy,l]
            elif not ix==jx and  iy==jy:
                for k in range(np.size(EI,1)):
                    Ke[i,j]+= mu*wi[k]*wj[jy]*l1df[ix,k]*l1df[jx,k]
    # 组装刚度矩阵
    for ele in EI:
        K[np.ix_(ele.flatten(), ele.flatten())]+=Ke
    print('K.size={}x{}'.format(numN,numN))
    return K
    
mass_inv = creat_mass_matrix(NC=NC,EI=EI)
K = creat_stiffness_matrix(NC=NC,EI=EI)
numN = np.size(NC,0)
# SE Solution, Time extrapolation
# ---------------------------------------------------------------
# initialize source time function and force vector f
# initialize time axis
t = np.arange(0, nt)*dt

# Initialization of the source time function
pt = 60*dt    # Gaussian width
t0 = 3*pt      # Time shift
src = -2/pt**2 * (t-t0) * np.exp(-1/pt**2 * (t-t0)**2)
# src  = ricker(dt,Tdom)
isrc = int(np.floor(numN/2))   # Source location
# isrc +=int( N/2+(nx*N+1)*(N/2))
# Initialization of solution vectors
u = np.zeros(numN, dtype =  float)
uold = np.zeros(numN, dtype =  float)
unew = np.zeros(numN, dtype =  float)
f = np.zeros(numN, dtype =  float) 
f[isrc]=1
print(u.shape)
px=[]
py=[]
for i in range(N*nx+1):
    px.append(NC[i,0])
for i in range(N*ny+1):
    py.append(NC[i*(N*nx+1),1])
xx,yy = np.meshgrid(px,py)
RemoveDir(filepath)
plt.figure()
for it in range(nt): 
    unew = dt**2 * mass_inv @ (f*src[it] - K @ u) + 2 * u - uold
    uold = u
    u = unew
    if not it % iplot:
        plt.pcolor(xx,yy,np.reshape(u,[N*nx+1,N*ny+1]))
        cb = plt.colorbar() 
        plt.axis('equal')
        plt.title('t={:.4f}s[{}x{}] N={} it={} vs={}'.format(dt*it,nx,ny,N,it,vs ))
        plt.savefig(filepath+'/{}.png'.format(it))
        cb.remove()
        print('{}%{}'.format(it, nt))


plt.show()
    

#!/usr/bin/env python
# -*- coding: utf-8 -*-
作者：@fm
import numpy as np


def legendre(N, x):
    """
    Returns the value of Legendre Polynomial P_N(x) at position x[-1, 1].
    """
    P = np.zeros(2 * N)

    if N == 0:
        P[0] = 1
    elif N == 1:
        P[1] = x
    else:
        P[0] = 1
        P[1] = x
    for i in range(2, N + 1):
        P[i] = (1.0 / float(i)) * ((2 * i - 1) * x * P[i - 1] - (i - 1) *
                                   P[i - 2])

    return(P[N])