Open phenomLi opened 5 years ago
在拥有了四叉树之后,我们就可以利用四叉树这个工具计算物体的N体力了,不过在计算N体力之前,我们首先要完成计算连接物体间的作用力。
连接物体间受到的力是胡克弹力,用作模拟某种弹性杆连接。胡克弹力的公式为: 其中k为弹簧系数,x0为弹簧原长,x为弹簧发生形变后的长度。
在对相互连接的物体添加弹力之前,我们首先要对ItemNode进行一些改造,使其支持物体连接:
ItemNode
// 最小区块节点(即叶节点,物体节点) export class ItemNode { //......省略 //连接的物体列表 public link: ItemNode[]; constructor(item) { //......省略 this.link = []; } /** * 计算弹力 * @param k 弹簧系数 */ calcElasticity(k: number) {} //......省略 }
代码中我们加入了一个link数组,用作存储与该物体连接的物体,然后,我们要在四叉树的主类QuadTree加入一个linkNode函数,用作连接两个物体:
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数组,对每个元素计算该元素产生的弹力即可:
calcElasticity
// 最小区块节点(即叶节点,物体节点) 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:
linkLength
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函数:
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所受到的合力,我们需要从四叉树的根节点开始,递归地执行如下步骤:
如果当前节点是一个叶子节点(而且它不是物体A),计算当前节点施加在物体A上的力,并将其加到A的合力上。
否则,计算商s/d的值。如果s/d<θ,将这个非叶子节点看成一个单独的物体,计算其施加在物体A上的力,并将其加到A的合力上。
否则,在当前节点的每个子节点上递归地执行上述步骤。
为了计算库伦力产生的N体力,我们在ItemNode类中增加一个calcElectromagnetic和coulombLaw方法用作计算库伦力:
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的源码,或者参考这篇文章。
在拥有了四叉树之后,我们就可以利用四叉树这个工具计算物体的N体力了,不过在计算N体力之前,我们首先要完成计算连接物体间的作用力。
胡克弹力
连接物体间受到的力是胡克弹力,用作模拟某种弹性杆连接。胡克弹力的公式为:
其中k为弹簧系数,x0为弹簧原长,x为弹簧发生形变后的长度。
在对相互连接的物体添加弹力之前,我们首先要对
ItemNode
进行一些改造,使其支持物体连接:代码中我们加入了一个
link
数组,用作存储与该物体连接的物体,然后,我们要在四叉树的主类QuadTree
加入一个linkNode
函数,用作连接两个物体:完成以上工作后,现在,我们就要编写
calcElasticity
函数的函数体了。思路很简单,遍历每个物体的link数组,对每个元素计算该元素产生的弹力即可:至于
linkLength
,我们暂时设定为30:库伦力
N体力在力导向图中体现为带电粒子的库伦力,库伦定律公式如下:
其中k为静电力常量,q1,q2为两个带电粒子的电荷量,r为粒子间距离,e方向向量。
接下来我们利用已有的四叉树进行物体间的库伦力计算,由于库伦力计算需要计算粒子间距离和粒子电荷量,所以首先我们给四叉树的两种节点加上 电荷量(charge) 属性:
叶子节点,也就是物体节点的电荷量默认为1,而非叶子节点的电荷量默认为0。
之后,我们就可以完善非叶子节点的
calcCentroid
函数和calcMass
函数:calcCentroid
函数会求当前非叶子节点代表的组的子节点的加权平均质心,calcMass
函数会求当前非叶子节点代表的组的质量和和电荷量和。接下来,重点来了,我们现在有了每个节点的质心和电荷量,可以使用四叉树求N体力了。
四叉树求N体力的核心是求非叶子节点和某叶子节点(物体)的受力关系,如果某个非叶子节点离某个物体并不足够远,那么就递归地遍历其所有子树。为了确定一个节点是否离得足够远,我们需要计算商 s/d, 其中s是非叶子节点所代表的区域的宽度,d是物体到节点所代表的那一组物体的质心的距离。然后将这个比值同一个阈值θ来作比较. 如果s/d<θ,那么非叶子节点就足够远。通过调整参数 θ,我们就可以来改变模拟的精度和速度。通常在实际中θ = 0.5是一个常常被使用的值。s和d关系如下图:![](https://github.com/phenomLi/Blog/raw/master/photos/2019-6-9/20161104233304006.png)
假如有一个为了物体A,那么为了计算物体A所受到的合力,我们需要从四叉树的根节点开始,递归地执行如下步骤:
如果当前节点是一个叶子节点(而且它不是物体A),计算当前节点施加在物体A上的力,并将其加到A的合力上。
否则,计算商s/d的值。如果s/d<θ,将这个非叶子节点看成一个单独的物体,计算其施加在物体A上的力,并将其加到A的合力上。
否则,在当前节点的每个子节点上递归地执行上述步骤。
为了计算库伦力产生的N体力,我们在
ItemNode
类中增加一个calcElectromagnetic
和coulombLaw
方法用作计算库伦力:库伦力的计算介绍到此即可。
最后
这两篇issue只围绕关于N体力的计算思路和实现,不会介绍如何实现力导向图,关于力导向图的绘制,大家可以阅读D3.js的源码,或者参考这篇文章。