N体受力问题（二）：计算物体作用力

[phenomLi]

在拥有了四叉树之后，我们就可以利用四叉树这个工具计算物体的N体力了，不过在计算N体力之前，我们首先要完成计算连接物体间的作用力。

胡克弹力

连接物体间受到的力是胡克弹力，用作模拟某种弹性杆连接。胡克弹力的公式为：

其中k为弹簧系数，x0为弹簧原长，x为弹簧发生形变后的长度。

在对相互连接的物体添加弹力之前，我们首先要对ItemNode进行一些改造，使其支持物体连接：

// 最小区块节点（即叶节点，物体节点）
export class ItemNode {

    //......省略

    //连接的物体列表
    public link: ItemNode[];

    constructor(item) {
        //......省略

        this.link = [];
    }

    /**
     * 计算弹力
     * @param k 弹簧系数
     */
    calcElasticity(k: number) {}
    
    //......省略
}

代码中我们加入了一个link数组，用作存储与该物体连接的物体，然后，我们要在四叉树的主类QuadTree加入一个linkNode函数，用作连接两个物体：

// 四叉树
export default class QuadTree {
    
    //...省略

    /**
     * 连接两个物体
     * @param node1 物体1
     * @param node2 物体2
     */
    linkNode(node1: ItemNode, node2: ItemNode) {
        // 相互加入对方的link数组
        node1.link.push(node2);
        node2.link.push(node1);
    }
}

完成以上工作后，现在，我们就要编写calcElasticity函数的函数体了。思路很简单，遍历每个物体的link数组，对每个元素计算该元素产生的弹力即可：

// 最小区块节点（即叶节点，物体节点）
export class ItemNode {

    //......省略

    /**
     * 计算弹力
     * @param k 弹簧系数
     */
    calcElasticity(k: number) {
        // 遍历link里的元素
        this.link.map(node => {
                // 获取该物体与连接物体间的质心连线向量
            let v = Vector.sub(this.centroid, node.centroid),
                // 计算弹力大小，其中linkLength表示连接杆的原始长度
                f = k*(Vector.len(v) - linkLength),
                // 弹力大小乘上质心连接向量，算出作用力
                ef = Vector.scl(f, Vector.nol(v));

            // 将计算出的力加到该物体的合外力
            this.force = Vector.add(this.force, ef);
        });
    }
    
    //......省略
}

至于linkLength，我们暂时设定为30：

export const linkLength = 30;

库伦力

N体力在力导向图中体现为带电粒子的库伦力，库伦定律公式如下：

其中k为静电力常量，q1，q2为两个带电粒子的电荷量，r为粒子间距离，e方向向量。

接下来我们利用已有的四叉树进行物体间的库伦力计算，由于库伦力计算需要计算粒子间距离和粒子电荷量，所以首先我们给四叉树的两种节点加上 电荷量（charge） 属性：

// 非最小区块节点
export class BlockNode {
    
    //......省略

    // 电荷量
    public charge: number;
    
    constructor(domain: Domain) {
        
        //......省略

        this.charge = 0;
    }

    //......省略
}



// 最小区块节点（即叶节点，物体节点）
export class ItemNode {
    //......省略

    // 电荷量
    public charge: number;
 
    constructor(item) {
        //......省略

        this.charge = 1;
    }
}

叶子节点，也就是物体节点的电荷量默认为1，而非叶子节点的电荷量默认为0。

之后，我们就可以完善非叶子节点的calcCentroid函数和calcMass函数：

// 非最小区块节点
export class BlockNode {
    
    //......省略

    // 该组质点计算
    calcCentroid() {
        let sumX = 0, sumY = 0;

        // 遍历子区块
        for(let block in this.childBlock) {
            let b = this.childBlock[block];

            // 若该子区块不为null
            if(b) {
                sumX += b.mass*b.centroid[0];
                sumY += b.mass*b.centroid[1];
            }
        }

        // 计算得质量加权质心
        this.centroid[0] = sumX/this.mass;
        this.centroid[1] = sumY/this.mass;
    }

    // 该组质量和电荷计算
    calcMass() {
        // 遍历子区块
        for(let block in this.childBlock) {
            // 若该子区块不为null
            if(this.childBlock[block]) {
                // 该节点的质量为子节点质量的和
                this.mass += this.childBlock[block].mass;
                // 该节点电荷量为子节点电荷量的和
                this.charge += this.childBlock[block].charge;
            }
        }
    }

    //......省略
}

calcCentroid函数会求当前非叶子节点代表的组的子节点的加权平均质心，calcMass函数会求当前非叶子节点代表的组的质量和和电荷量和。

接下来，重点来了，我们现在有了每个节点的质心和电荷量，可以使用四叉树求N体力了。

四叉树求N体力的核心是求非叶子节点和某叶子节点（物体）的受力关系，如果某个非叶子节点离某个物体并不足够远，那么就递归地遍历其所有子树。为了确定一个节点是否离得足够远，我们需要计算商 s/d, 其中s是非叶子节点所代表的区域的宽度，d是物体到节点所代表的那一组物体的质心的距离。然后将这个比值同一个阈值θ来作比较. 如果s/d<θ，那么非叶子节点就足够远。通过调整参数 θ，我们就可以来改变模拟的精度和速度。通常在实际中θ = 0.5是一个常常被使用的值。s和d关系如下图：

假如有一个为了物体A，那么为了计算物体A所受到的合力，我们需要从四叉树的根节点开始，递归地执行如下步骤：

1. 如果当前节点是一个叶子节点（而且它不是物体A），计算当前节点施加在物体A上的力，并将其加到A的合力上。

2. 否则，计算商s/d的值。如果s/d<θ，将这个非叶子节点看成一个单独的物体，计算其施加在物体A上的力，并将其加到A的合力上。

3. 否则，在当前节点的每个子节点上递归地执行上述步骤。

为了计算库伦力产生的N体力，我们在ItemNode类中增加一个calcElectromagnetic和coulombLaw方法用作计算库伦力：

// 最小区块节点（即叶节点，物体节点）
export class ItemNode {

    //......省略

    /**
     * 计算电磁力（库伦力）
     * @param root 当前节点
     * @param k 静电力常量
     */
    calcElectromagnetic(root: NodeType, k: number) {
        // θ
        let θ = 0.5;

        // 若是非叶子节点
        if(root instanceof BlockNode) {
            // 计算 s/d
            let n = root.domain.width/Vector.len(Vector.sub(root.centroid, this.centroid));

            // s/d < θ：足够远，看成一组物体
            if(n < θ) {
                // 套用库伦定律公式计算该节点和该组间的库伦力，并将结果加到该节点的合力上
                this.force = Vector.add(this.force, this.coulombLaw(this, root, k)); 
            }
            // s/d >= θ：不够远，不看成一组物体，则继续递归计算
            else {
                // 遍历子节点
                for(let block in root.childBlock) {
                    // 若子节点不为null
                    if(root.childBlock[block]) {
                        // 递归计算
                        this.calcElectromagnetic(root.childBlock[block], k);
                    }
                }
            }

        }

        // 若是叶子节点
        if(root instanceof ItemNode) {
            // 直接套用库伦定律公式计算库伦力，并将结果加到该节点的合力上
            this.force = Vector.add(this.force, this.coulombLaw(this, root, k)); 
        }
    }

    /**
     * 库伦定律计算公式
     * @param node1 节点1
     * @param node2 节点2
     * @param k 静电力常量
     */
    coulombLaw(node1: ItemNode, node2: NodeType, k: number): vector {
        let v = Vector.sub(node1.centroid, node2.centroid),
                len = Vector.len(v);

        return Vector.scl(k*node1.charge*node2.charge/(len*len), Vector.nol(v));
    }
    
    //......省略
}

库伦力的计算介绍到此即可。

最后

这两篇issue只围绕关于N体力的计算思路和实现，不会介绍如何实现力导向图，关于力导向图的绘制，大家可以阅读D3.js的源码，或者参考这篇文章。