This module defines the trait FiniteElement, and provides methods to solve finite element systems.

Example

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);

Use of the macro provided by finiteelement_macro

To use one of the macro provided by finiteelement_marco, call it without parameter at the begining or your code. The piece of code surrounded by //======== in the example above can be replaced by a call to 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);
}