2 个不稳定版本
0.3.0 | 2019年8月29日 |
---|---|
0.1.0 | 2019年7月15日 |
#11 in #norm
210KB
4.5K SLoC
此模块定义了 FiniteElement
特征,并提供了解决有限元系统的方法。
示例
use finiteelement::*;
//Define a struct that implements the FiniteElement trait
#[derive(Clone)]
pub struct Spring {
pub a: usize,
pub b: usize,
pub l: f64,
pub k: f64,
}
impl FiniteElement<f64> for Spring {
fn forces(&self, positions: &[Point<f64>], forces: &mut [f64]) {
// add to both a and b the force resulting from this spring.
let ab = positions[self.b].clone() - positions[self.a].clone();
let norm = ab.norm();
forces[3 * self.a] += self.k * (norm - self.l) * ab.x / norm;
forces[3 * self.a + 1] += self.k * (norm - self.l) * ab.y / norm;
forces[3 * self.a + 2] += self.k * (norm - self.l) * ab.z / norm;
forces[3 * self.b] -= self.k * (norm - self.l) * ab.x / norm;
forces[3 * self.b + 1] -= self.k * (norm - self.l) * ab.y / norm;
forces[3 * self.b + 2] -= self.k * (norm - self.l) * ab.z / norm;
}
fn jacobian(&self, positions: &[Point<f64>], jacobian: &mut Sparse<f64>) {
// add to both a and b the force resulting from this self.
let ab = positions[self.b].clone() - positions[self.a].clone();
let norm = ab.norm();
let norm3 = norm * norm * norm;
for u in 0..3 {
for v in 0..3 {
let j = if u == v {
self.k * (1. - self.l / norm + self.l * ab[u] * ab[u] / norm3)
} else {
self.k * self.l * ab[u] * ab[v] / norm3
};
// Change in the force on B when moving B.
jacobian[3 * self.b + u][3 * self.b + v] -= j;
// Change in the force on A when moving B.
jacobian[3 * self.a + u][3 * self.b + v] += j;
// Change in the force on B when moving A.
jacobian[3 * self.b + u][3 * self.a + v] += j;
// Change in the force on A when moving A.
jacobian[3 * self.a + u][3 * self.a + v] -= j;
}
}
}
}
let elts = [
Spring {
a: 0,
b: 1,
l: 1.,
k: 1.,
},
Spring {
a: 1,
b: 2,
l: 2.,
k: 0.5,
},
Spring {
a: 0,
b: 2,
l: 3.,
k: 5.
},
];
let system = (0..elts.len()).map(|i| {Box::new(elts[i].clone()) as Box<dyn FiniteElement<f64>>}).collect::<Vec<_>>();
let mut positions = vec![
Point { x: 0., y: 0., z: 0. },
Point {x : 1., y: 0., z: 0.},
Point{x: 0., y: 1., z: 1.}];
let mut ws = FesWorkspace::new(positions.len());
let epsilon_stop = 1e-4;
let gradient_switch = 1e-3;
let mut solved = false;
for i in (0..20) {
solved = fes_one_step(&system, &mut positions, epsilon_stop, gradient_switch, &mut ws);
if solved {
break;
}
}
assert!(solved);
finiteelement_macro 提供的宏的使用
要使用 finiteelement_marco
提供的宏之一,请在代码开头不带参数调用它。上面示例中用 //========
包围的代码可以替换为对 finiteelement_macros::auto_impl_spring!{}
的调用
use finiteelement::*;
use std::borrow::Borrow;
auto_impl_spring!{}
pub fn main() {
let elts = [
Spring {
a: 0,
b: 1,
l: 1.,
k: 1.,
},
Spring {
a: 1,
b: 2,
l: 2.,
k: 0.5,
},
Spring {
a: 0,
b: 2,
l: 3.,
k: 5.
},
];
let system = (0..elts.len()).map(|i| {Box::new(elts[i].clone()) as Box<dyn FiniteElement<f64>>}).collect::<Vec<_>>();
let mut positions = vec![
Point { x: 0., y: 0., z: 0. },
Point {x : 1., y: 0., z: 0.},
Point{x: 0., y: 1., z: 1.}];
let mut ws = FesWorkspace::new(positions.len());
let epsilon_stop = 1e-4;
let gradient_switch = 1e-3;
let mut solved = false;
for i in (0..20) {
solved = fes_one_step(&system, &mut positions, epsilon_stop, gradient_switch, &mut ws);
if solved {
break;
}
}
assert!(solved);
}
依赖
~1.5–2.5MB
~49K SLoC