RBF三维应力插值

算法原理

参见上一篇：径向基函数（RBF）插值

测试代码

下文简单测试了三维下的RBF函数插值。

测试代码如下：

Test_RBF_interp3d.py

import numpy as np
import pandas as pd
from sklearn.neighbors import KNeighborsClassifier  # KNN
from random import choice
from scipy.interpolate import Rbf

import RBF
import visualization

if __name__ == "__main__":
    data = pd.read_csv("ES_result.csv").values
    all_points, all_stress = data[-1000:, :3], data[-1000:, 3]  # 共计1000个节点
    knn = KNeighborsClassifier(n_neighbors = 1, algorithm = 'ball_tree')
    knn.fit(all_points, [1] * len(all_points))  # 构建KNN
    indexes = list(range(len(all_points)))
    sample_points = []
    distance, n_points = knn.kneighbors(all_points, n_neighbors = 2, return_distance = True)  # 寻找邻近的15个点
    while len(indexes) >= 15:
        sample_index = choice(indexes)
        sample_points.append(np.mean(all_points[[n_points[sample_index]], :].reshape(-1, 3), axis = 0))  # 按列求平均值构造采样点
        indexes = [i for i in indexes if i not in n_points[sample_index]]
    sample_points = np.array(sample_points)
    # print(sample_points.shape)
    # visualization.viz_nodes(sample_points)
    w = RBF.rbf_coefficient(all_points, all_stress, "linear")
    pre_stress = RBF.rbf_interpolation(all_points, w, sample_points, "linear")
    # funs = Rbf(all_points[:, 0], all_points[:, 1], all_points[:, 2], all_stress, function = 'gaussian')
    # pre_stress = funs(sample_points[:, 0], sample_points[:, 1], sample_points[:, 2])
    # print(pre_stress)
    visualization.viz_nodes(np.vstack((all_points, sample_points)),
                            np.vstack((all_stress.reshape(-1, 1), pre_stress.reshape(-1, 1))))
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31

测试结果

结果中，小点点为控制节点坐标，小叉叉为插值节点坐标。节点颜色的深浅与应力大小相映射。

• gaussian基函数

  

• cubic基函数

  

• linear基函数

  

总结

• 高斯核宽度取值影响插值结果

  当基函数为高斯基函数时，其核函数如下：
  $$\phi (r)=\exp \left(-\frac{r^2}{2{\sigma }^2}\right)$$
  其中 $r$ 为欧式距离，$\sigma$ 为核函数宽度。

  根据测试可知，$\sigma$ 取值不同将得到不同的插值结果。鉴于scipy库中已实现RBF插值函数的封装，打开其函数定义，我们可以发现，在scipy中，高斯基函数定义如下：

  def _h_gaussian(self, r):
      return np.exp(-(1.0/self.epsilon*r)**2)
  1
  2

  即
  $$\phi (r)=\exp \left(-\frac{r^2}{{ϵ}^2}\right)$$
  scipy插值库中self.epsilon默认取基于边界超立方体的“节点之间的平均距离”，计算代码如下：

  if self.epsilon is None:
  # default epsilon is the "the average distance between nodes" based on a bounding hypercube
  ximax = np.amax(self.xi, axis=1)
  ximin = np.amin(self.xi, axis=1)
  edges = ximax - ximin
  edges = edges[np.nonzero(edges)]
  self.epsilon = np.power(np.prod(edges)/self.N, 1.0/edges.size)
  1
  2
  3
  4
  5
  6
  7

  在本三维插值测试中，$ϵ$ 的计算公式为
  $$ϵ=\sqrt[3]{\frac{V}{N}}$$
  其中 $V$ 为所有节点包围盒的体积，$N$ 为节点总数。

  仿照scipy源码，我们可以如下计算得到高斯核宽度
  $$\sigma =\frac{\sqrt{2}}{2}\sqrt[3]{\frac{V}{N}}$$

  对于Multiquadrics基函数中的 $\sigma$ 也应当合理取值。