（二）、 Fealpy 重心坐标导数的计算

我们先给出三角形单元的中心坐标公式

对应的代码为：

import numpy as np
from fealpy.mesh import MeshFactory
mf = MeshFactory()
box = [0, 1, 0, 1]
mesh = mf.boxmesh2d(box, nx=1, ny=1, meshtype='tri')
NC = mesh.number_of_cells()

node = mesh.entity('node')
cell = mesh.entity('cell')
v0 = node[cell[:, 2], :] - node[cell[:, 1], :]  # x2 - x1 （NC， 2）
v1 = node[cell[:, 0], :] - node[cell[:, 2], :]  # x0 - x2
v2 = node[cell[:, 1], :] - node[cell[:, 0], :]  # x1 - x0
nv = np.cross(v2, -v1)  # 叉乘  2倍三角形面积

Dlambda = np.zeros((NC, 3, 2), dtype=np.float64)  # 单元 * 重心坐标数量 * 每个坐标分量
length= nv  # (NC,)
W = np.array([[0, 1], [-1, 0]], dtype=np.int_)
Dlambda[:, 0, :] = v0@W/length.reshape(-1, 1)
Dlambda[:, 1, :] = v1@W/length.reshape(-1, 1)
Dlambda[:, 2, :] = v2@W/length.reshape(-1, 1)

计算一维区间的导数

对应的代码为：

import numpy as np
from fealpy.mesh import IntervalMesh
node = np.array([[0.0], [0.5], [1.0]], dtype=np.float64)
cell = np.array([[0, 1], [1, 2]], dtype=np.int_)
mesh = IntervalMesh(node, cell)
NC = mesh.number_of_cells()
v0 = node[cell[:, 1]] - node[cell[:, 0]]  # x1 - x0
Dlambda = np.zeros((NC, 2, 1), dtype=np.float64)
Dlambda[:, 0, :] = 1 / v0
Dlambda[:, 1, :] = -1 / v0