# 重心坐标及其导数

$$(\lambda_1+\cdots+\lambda_n)p=\lambda_1v_1+\cdots+\lambda_nv_n$$

## 一维重心坐标

### 一维重心坐标求导

$$\left[\begin{matrix}x_0 & x_1\\1 & 1\end{matrix}\right]\left[\begin{matrix}\lambda_0\\ \lambda_1\end{matrix}\right]=\left[\begin{matrix}x\\1\end{matrix}\right]$$

$$\lambda_0=\frac{x_1-x}{x_1-x_0},\quad \lambda_1=\frac{x-x_0}{x_1-x_0}$$

$$\frac{d\lambda_0}{dx}=\frac{-1}{x_1-x_0},\quad \frac{d\lambda_1}{dx}=\frac{1}{x_1-x_0}$$

### 一维重心坐标的物理意义

import numpy as np
from fealpy.mesh import IntervalMesh
import matplotlib.pyplot as plt

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)
fig = plt.figure()
axes = fig.gca()
mesh.find_node(axes, showindex = True)
mesh.find_cell(axes, showindex = True)
plt.show()

number_of_cells = mesh.number_of_cells()
node = mesh.entity('node')
cell = mesh.entity('cell')

length_of_cells = node[cell[:, 1], :] - node[cell[:, 0], :]
Dlambda = np.zeros((number_of_cells, 2, 1), dtype = np.float64)

Dlambda[:, 0, :] = -1 / length_of_cells
Dlambda[:, 1, :] = 1 / length_of_cells
print(Dlambda)

## 二维重心坐标

### 二维重心坐标函数求导

$$A\lambda=X$$

$$\left[\begin{matrix}x_{00} & x_{10}& x_{20}\\x_{01}&x_{11}&x_{21}\\1&1 & 1\end{matrix}\right]\left[\begin{matrix}\lambda_0\\ \lambda_1\\\lambda_2\end{matrix}\right]=\left[\begin{matrix}x_0\\x_1\\1\end{matrix}\right]$$

$$A^{-1}=\left[\begin{matrix}a_{00}&a_{01}&a_{02}\\a_{10}&a_{11}&a_{12}\\a_{20}&a_{21}&a_{22}\end{matrix}\right]$$

$$\left[\begin{matrix}\lambda_0\\ \lambda_1\\\lambda_2\end{matrix}\right]=\left[\begin{matrix}a_{00}&a_{01}&a_{02}\\a_{10}&a_{11}&a_{12}\\a_{20}&a_{21}&a_{22}\end{matrix}\right]\left[\begin{matrix}x_0\\x_1\\1\end{matrix}\right]$$

$$\nabla\lambda_0=\left[\begin{matrix}\frac{\partial\lambda_0}{\partial x_0}\\\frac{\partial\lambda_0}{\partial x_1}\\\end{matrix}\right]=\left[\begin{matrix}a_{00}\\a_{01}\end{matrix}\right]$$

$$\nabla\lambda_1=\left[\begin{matrix}\frac{\partial\lambda_1}{\partial x_0}\\\frac{\partial\lambda_1}{\partial x_1}\\\end{matrix}\right]=\left[\begin{matrix}a_{10}\\a_{11}\end{matrix}\right]$$

$$\nabla\lambda_2=\left[\begin{matrix}\frac{\partial\lambda_2}{\partial x_0}\\\frac{\partial\lambda_2}{\partial x_1}\\\end{matrix}\right]=\left[\begin{matrix}a_{20}\\a_{21}\end{matrix}\right]$$

### 二维重心坐标的物理意义

• $\lambda_0$表示$x_0$所对的三角形$\triangle_{x,x_1,x_2}$占整个三角形$\triangle_{x_0,x_1,x_2}$面积的比例。
• $\lambda_1$表示$x_1$所对的三角形$\triangle_{x,x_2,x_0}$占整个三角形$\triangle_{x_0,x_1,x_2}$面积的比例。
• $\lambda_0$表示$x_2$所对的三角形$\triangle_{x,x_0,x_1}$占整个三角形$\triangle_{x_0,x_1,x_2}$面积的比例。

import numpy as np
from fealpy.mesh import MeshFactory as MF

box = [0, 1, 0, 1]
mesh = MF.boxmesh2d(box, nx = 1, ny = 1, meshtype = 'tri')
number_of_cells = mesh.number_of_cells()

import matplotlib.pyplot as plt
fig = plt.figure()
axes = fig.gca()
mesh.find_node(axes, showindex = True)
mesh.find_cell(axes, showindex = True)
plt.show()

node = mesh.entity('node')
cells = mesh.entity('cell')
# print(cells.shape)

i = 0
Dlambda = np.zeros((number_of_cells, 3, 2), dtype=float)
for cell in cells:
# print(cell)
transform_matrix_A = np.array([node[cell, 0], node[cell, 1], np.ones((3), dtype = np.float64)], dtype = np.float64)
# print(transform_matrix_A)
inv_A = np.linalg.inv(transform_matrix_A)
Dlambda[i, 0, :] = np.array([inv_A[0, 0], inv_A[0, 1]], dtype = float)
Dlambda[i, 1, :] = np.array([inv_A[1, 0], inv_A[1, 1]], dtype = float)
Dlambda[i, 2, :] = np.array([inv_A[2, 0], inv_A[2, 1]], dtype = float)
i += 1
# print(i, " cell:\n", transform_matrix_A)
# print(inv_A)
print(Dlambda)

### 向量化编程

import numpy as np
from fealpy.mesh import MeshFactory as MF
import matplotlib.pyplot as plt

box = [0, 1, 0, 1]
mesh = MF.boxmesh2d(box, nx = 1, ny = 1, meshtype = 'tri')

number_of_cells = mesh.number_of_cells()
node = mesh.entity('node')
cells = mesh.entity('cell')

cells_A = np.zeros((3, 3, number_of_cells), dtype = float)
inv_cells_A = np.zeros((3, 3, number_of_cells), dtype = float)
cells_A[:, :, 0: number_of_cells] = np.array([node[cells[0: number_of_cells], 0].T, node[cells[0: number_of_cells], 1].T, np.ones_like(node[cells[0: number_of_cells], 0], dtype = float).T], dtype = float)
inv_cells_A = np.array([np.linalg.inv(cells_A[:, :, i]) for i in range(number_of_cells)], dtype = float)
Dlambda = inv_cells_A[:, :, 0: 2]
print(Dlambda)

# NC = mesh.number_of_cells()
# node = mesh.entity('node')
# cell = mesh.entity('cell')
# v0 = node[cell[:, 2], :] - node[cell[:, 1], :] # x 2 − x 1
# v1 = node[cell[:, 0], :] - node[cell[:, 2], :] # x 0 − x 2
# v2 = node[cell[:, 1], :] - node[cell[:, 0], :] # x 1 − x 0
# nv = np.cross(v2, -v1)

# Dlambda = np.zeros((NC, 3, 2), dtype=np.float64)
# length = nv #
# 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)
# print(Dlambda)

## 三维重心坐标函数

### 三维重心坐标函数求导

$$A\lambda=X$$

$$\left[\begin{matrix}x_{00} & x_{10}& x_{20}& x_{30}\\x_{01}&x_{11}&x_{21}& x_{31}\\x_{02}&x_{12}&x_{22}&x_{32}\\1&1 & 1&1\end{matrix}\right]\left[\begin{matrix}\lambda_0\\ \lambda_1\\\lambda_2\\\lambda_3\end{matrix}\right]=\left[\begin{matrix}x_0\\x_1\\x_2\\1\end{matrix}\right]$$

$$\lambda=A^{-1}x$$

$$A^{-1}=\left[\begin{matrix}a_{00}&a_{01}&a_{02}&a_{03}\\a_{10}&a_{11}&a_{12}&a_{13}\\a_{20}&a_{21}&a_{22}&a_{23}\\a_{30}&a_{31}&a_{32}&a_{33}\end{matrix}\right]$$

$$\left[\begin{matrix}\lambda_0\\ \lambda_1\\\lambda_2\\\lambda_3\end{matrix}\right]=\left[\begin{matrix}a_{00}&a_{01}&a_{02}&a_{03}\\a_{10}&a_{11}&a_{12}&a_{13}\\a_{20}&a_{21}&a_{22}&a_{23}\\a_{30}&a_{31}&a_{32}&a_{33}\end{matrix}\right]\left[\begin{matrix}x_0\\x_1\\x_2\\1\end{matrix}\right]$$

$$\nabla\lambda_0=\left[\begin{matrix}\frac{\partial\lambda_0}{\partial x_0}\\\frac{\partial\lambda_0}{\partial x_1}\\\frac{\partial\lambda_0}{\partial x_2}\\\end{matrix}\right]=\left[\begin{matrix}a_{00}\\a_{01}\\a_{02}\end{matrix}\right]$$

$$\nabla\lambda_1=\left[\begin{matrix}\frac{\partial\lambda_1}{\partial x_0}\\\frac{\partial\lambda_1}{\partial x_1}\\\frac{\partial\lambda_1}{\partial x_2}\\\end{matrix}\right]=\left[\begin{matrix}a_{10}\\a_{11}\\a_{12}\end{matrix}\right]$$

$$\nabla\lambda_2=\left[\begin{matrix}\frac{\partial\lambda_2}{\partial x_0}\\\frac{\partial\lambda_2}{\partial x_1}\\\frac{\partial\lambda_2}{\partial x_2}\\\end{matrix}\right]=\left[\begin{matrix}a_{20}\\a_{21}\\a_{22}\end{matrix}\right]$$

$$\nabla\lambda_3=\left[\begin{matrix}\frac{\partial\lambda_3}{\partial x_0}\\\frac{\partial\lambda_3}{\partial x_1}\\\frac{\partial\lambda_3}{\partial x_2}\\\end{matrix}\right]=\left[\begin{matrix}a_{30}\\a_{31}\\a_{32}\end{matrix}\right]$$

### 三维重心坐标的物理意义

• $\lambda_0$表示$x_0$所对的四面体$\triangle_{x,x_1,x_2,x_3}$占整个四面体$\triangle_{x_0,x_1,x_2，x_3}$体积的比例。
• $\lambda_1$表示$x_1$所对的四面体$\triangle_{x,x_0,x_2,x_3}$占整个四面体$\triangle_{x_0,x_1,x_2，x_3}$体积的比例。
• $\lambda_2$表示$x_2$所对的四面体$\triangle_{x,x_0,x_1,x_3}$占整个四面体$\triangle_{x_0,x_1,x_2，x_3}$体积的比例。
• $\lambda_3$表示$x_3$所对的四面体$\triangle_{x,x_0,x_1,x_2}$占整个四面体$\triangle_{x_0,x_1,x_2，x_3}$体积的比例。

import numpy as np
from fealpy.mesh import MeshFactory as MF
import matplotlib.pyplot as plt

mesh = MF.one_tetrahedron_mesh(meshtype = 'iso')
number_of_cells = mesh.number_of_cells()

fig = plt.figure()
axes = fig.gca(projection='3d')
mesh.find_node(axes, showindex=True)
mesh.find_cell(axes, showindex=True)
plt.show()

node = mesh.entity('node')
cells = mesh.entity('cell')
print(cells.shape)

i = 0
Dlambda = np.zeros((number_of_cells, 4, 3), dtype=float)
for cell in cells:
print(cell)
transform_matrix_A = np.array([node[cell, 0], node[cell, 1], node[cell, 2], np.ones((4), dtype = np.float64)], dtype = np.float64)
inv_A = np.linalg.inv(transform_matrix_A)
Dlambda[i, 0, :] = np.array([inv_A[0, 0], inv_A[0, 1], inv_A[0, 2]], dtype = float)
Dlambda[i, 1, :] = np.array([inv_A[1, 0], inv_A[1, 1], inv_A[1, 2]], dtype = float)
Dlambda[i, 2, :] = np.array([inv_A[2, 0], inv_A[2, 1], inv_A[2, 2]], dtype = float)
Dlambda[i, 3, :] = np.array([inv_A[3, 0], inv_A[3, 1], inv_A[3, 2]], dtype = float)
i += 1
# print(i, " cell:\n", transform_matrix_A)
# print(inv_A)
print(Dlambda)

### 向量化编程

import numpy as np
from fealpy.mesh import MeshFactory as MF
import matplotlib.pyplot as plt

mesh = MF.one_tetrahedron_mesh(meshtype = 'iso')
number_of_cells = mesh.number_of_cells()

fig = plt.figure()
axes = fig.gca(projection='3d')
mesh.find_node(axes, showindex=True)
mesh.find_cell(axes, showindex=True)
plt.show()

node = mesh.entity('node')
cells = mesh.entity('cell')
print(cells.shape)

cells_A = np.zeros((4, 4, number_of_cells), dtype = float)
inv_cells_A = np.zeros((4, 4, number_of_cells), dtype = float)
cells_A[:, :, 0: number_of_cells] = np.array([node[cells[0: number_of_cells], 0].T, node[cells[0: number_of_cells], 1].T, node[cells[0: number_of_cells], 2].T, np.ones_like(node[cells[0: number_of_cells], 0], dtype = float).T], dtype = float)
inv_cells_A = np.array([np.linalg.inv(cells_A[:, :, i]) for i in range(number_of_cells)], dtype = float)
Dlambda = inv_cells_A[:, :, 0: 3]
print(Dlambda)