OpenTissue::fem::detail Namespace Reference

Classes
class	NodeTraits
class	TetrahedronTraits
Functions
template<typename tetrahedron_iterator , typename real_type >
void	add_plasticity_force (tetrahedron_iterator begin, tetrahedron_iterator end, real_type const &dt)
template<typename node_iterator >
void	clear_stiffness_assembly (node_iterator begin, node_iterator end)
template<typename real_type , typename vector3_type >
void	compute_B (vector3_type const &e10, vector3_type const &e20, vector3_type const &e30, real_type volume, vector3_type *B)
template<typename real_type , typename vector3_type >
void	compute_isotropic_elasticity_vector (real_type const &young, real_type const &poisson, vector3_type &D)
template<typename real_type , typename matrix_type >
void	compute_isotropic_elasticity_matix (real_type const &young, real_type const &poisson, matrix_type &D)
template<typename real_type >
void	compute_isotropic_elasticity_vector (real_type const &young, real_type const &poisson, real_type &D00, real_type &D01, real_type &D33)
template<typename real_type , typename vector3_type , typename matrix3x3_type >
void	compute_Ke (vector3_type const &p0, vector3_type const &p1, vector3_type const &p2, vector3_type const &p3, real_type const &E, real_type const &nu, matrix3x3_type Ke[4][4])
template<typename real_type , typename vector3_type , typename matrix3x3_type >
void	compute_Ke (vector3_type *B, vector3_type const &D, real_type const &volume, matrix3x3_type Ke[4][4])
template<typename real_type , typename vector3_type , typename matrix3x3_type >
void	compute_Ke_ij (vector3_type const &Bi, vector3_type const &D, vector3_type const &Bj, real_type const &volume, matrix3x3_type &Ke_ij)
template<typename fem_mesh >
void	compute_mass (fem_mesh &mesh)
template<typename vector3_type >
vector3_type::value_type	compute_volume (vector3_type const &e10, vector3_type const &e20, vector3_type const &e30)
template<typename vector3_type >
vector3_type::value_type	compute_volume (vector3_type const &p0, vector3_type const &p1, vector3_type const &p2, vector3_type const &p3)
template<typename fem_mesh >
void	conjugate_gradients (fem_mesh &mesh, unsigned int min_iterations, unsigned int max_iterations)
template<typename fem_mesh , typename real_type >
void	dynamics_assembly (fem_mesh &mesh, real_type const &mass_damping, real_type const &dt)
template<typename tetrahedron_iterator , typename real_type >
void	initialize_plastic (tetrahedron_iterator begin, tetrahedron_iterator end, real_type const &c_yield, real_type const &c_creep, real_type const &c_max)
template<typename tetrahedron_iterator >
void	initialize_stiffness_elements (tetrahedron_iterator begin, tetrahedron_iterator end)
template<typename fem_mesh , typename real_type >
void	position_update (fem_mesh &mesh, real_type const &dt)
template<typename tetrahedron_iterator >
void	reset_orientation (tetrahedron_iterator const &begin, tetrahedron_iterator const &end)
template<typename tetrahedron_iterator >
void	stiffness_assembly (tetrahedron_iterator const &begin, tetrahedron_iterator const &end)
template<typename real_type , typename tetrahedron_iterator >
void	uniform_density (tetrahedron_iterator const &begin, tetrahedron_iterator const &end, real_type const &density)
template<typename real_type , typename tetrahedron_iterator >
void	uniform_poisson (tetrahedron_iterator const &begin, tetrahedron_iterator const &end, real_type const &poisson)
template<typename real_type , typename tetrahedron_iterator >
void	uniform_young (tetrahedron_iterator const &begin, tetrahedron_iterator const &end, real_type const &young)
template<typename tetrahedron_iterator >
void	update_orientation (tetrahedron_iterator const &begin, tetrahedron_iterator const &end)

---

Function Documentation

template<typename tetrahedron_iterator , typename real_type >

void OpenTissue::fem::detail::add_plasticity_force	(	tetrahedron_iterator	begin,
		tetrahedron_iterator	end,
		real_type const &	dt
	)			[inline]

Add plasticity forces.

Parameters:

	begin	Iterator to first tetrahedron.
	end	Iterator to one past last tetrahedron.
	dt	Simulation time step

template<typename node_iterator >

void OpenTissue::fem::detail::clear_stiffness_assembly	(	node_iterator	begin,
		node_iterator	end
	)			[inline]

Parameters:

	begin
	end

template<typename real_type , typename vector3_type >

void OpenTissue::fem::detail::compute_B	(	vector3_type const &	e10,
		vector3_type const &	e20,
		vector3_type const &	e30,
		real_type	volume,
		vector3_type *	B
	)			[inline]

Calculate B-matrices (derivatives of shape functions N, i.e B = SN).

Parameters:

	e10	Given the four corner points p0,p1,p2, and p3. This is p1-p0.
	e20	Given the four corner points p0,p1,p2, and p3. This is p2-p0.
	e30	Given the four corner points p0,p1,p2, and p3. This is p3-p0.
	volume	The volume of the tetrahedron.
	B	Upon return this array of 4 B-vectors contains the computed values.

template<typename real_type , typename matrix_type >

void OpenTissue::fem::detail::compute_isotropic_elasticity_matix	(	real_type const &	young,
		real_type const &	poisson,
		matrix_type &	D
	)			[inline]

Compute Elasticy Matrix.

Parameters:

	young	Youngs modulus (stiffness)
	poisson	Poissons ratio (compressability)
	D	Upon return holds the isotropoc linear elasticity matrix.

template<typename real_type , typename vector3_type >

void OpenTissue::fem::detail::compute_isotropic_elasticity_vector	(	real_type const &	young,
		real_type const &	poisson,
		vector3_type &	D
	)			[inline]

Compute Isotropic Elasticity Matrix.

Parameters:

	young	Young Modulus.
	poisson	Poisson ratio.
	D	Upon return holds the elasticity matrix in vector form D = [D0,D1,D2]

template<typename real_type >

void OpenTissue::fem::detail::compute_isotropic_elasticity_vector	(	real_type const &	young,
		real_type const &	poisson,
		real_type &	D00,
		real_type &	D01,
		real_type &	D33
	)			[inline]

Alternative representation of isotropic elasticity matrix.

Parameters:

	young	Youngs modulus (stiffness)
	poisson	Poissons ratio (compressability)
	D00
	D01
	D33

template<typename real_type , typename vector3_type , typename matrix3x3_type >

void OpenTissue::fem::detail::compute_Ke	(	vector3_type *	B,
		vector3_type const &	D,
		real_type const &	volume,
		matrix3x3_type	Ke[4][4]
	)			[inline]

Compute Stiffness matrix of tetrahedral element.

template<typename real_type , typename vector3_type , typename matrix3x3_type >

void OpenTissue::fem::detail::compute_Ke	(	vector3_type const &	p0,
		vector3_type const &	p1,
		vector3_type const &	p2,
		vector3_type const &	p3,
		real_type const &	E,
		real_type const &	nu,
		matrix3x3_type	Ke[4][4]
	)			[inline]

Compute Stiffness matrix of tetrahedral element.

Parameters:

	p0
	p1
	p2
	p3
	E	Youngs modulus.
	nu	Poisson ratio.
	Ke	Upon return contains the computed stiffness values.

template<typename real_type , typename vector3_type , typename matrix3x3_type >

void OpenTissue::fem::detail::compute_Ke_ij	(	vector3_type const &	Bi,
		vector3_type const &	D,
		vector3_type const &	Bj,
		real_type const &	volume,
		matrix3x3_type &	Ke_ij
	)			[inline]

Compute i,j sub-block of the element stiffness matrix Ke

Parameters:

	Bi	Sub-block entries of B matrix of i'th node given in vector-form. Ie. Bi = [bi ci di]^T.
	D	The elasticity matrix in vector form D = [D0 D1 D2]^T.
	Bj	Sub-block entries of B matrix of i'th node given in vector-form. Ie. Bj = [bj cj dj]^T.
	V	The volume of the tetrahedron.
	Ke_ij	Upon return, contains the computed value.

template<typename fem_mesh >

void OpenTissue::fem::detail::compute_mass

(

fem_mesh &

mesh

)

[inline]

Compute Mass.

Note: Volume should have been computed prior to invoking this function. Thus make sure that initialize_stiffness_elements have been invoked prior to this function.

Parameters:

mesh

template<typename vector3_type >

vector3_type::value_type OpenTissue::fem::detail::compute_volume	(	vector3_type const &	e10,
		vector3_type const &	e20,
		vector3_type const &	e30
	)			[inline]

Compute volume of tetrahedron.

Parameters:

	e10	edge vector from p0 to p1, i.e. e10 = p1-p0.
	e20	edge vector from p0 to p2, i.e. e20 = p2-p0.
	e30	edge vector from p0 to p3, i.e. e30 = p3-p0.

Returns:

The signed volume of the tetrahedra with nodes p0,p1,p2 and p3.

template<typename vector3_type >

vector3_type::value_type OpenTissue::fem::detail::compute_volume	(	vector3_type const &	p0,
		vector3_type const &	p1,
		vector3_type const &	p2,
		vector3_type const &	p3
	)			[inline]

Compute Tetrahedron Volume

Parameters:

	p0
	p1
	p2
	p3

Returns:

The signed volume of the tetrahedra with nodes p0,p1,p2 and p3.

template<typename fem_mesh >

void OpenTissue::fem::detail::conjugate_gradients	(	fem_mesh &	mesh,
		unsigned int	min_iterations,
		unsigned int	max_iterations
	)			[inline]

Conjugate Gradient Solver. This solver has been hard-wired for a WarpT4Mesh to solve the equation

A v = b

template<typename fem_mesh , typename real_type >

void OpenTissue::fem::detail::dynamics_assembly	(	fem_mesh &	mesh,
		real_type const &	mass_damping,
		real_type const &	dt
	)			[inline]

Setup dynamic equation. This is the equation of motion given as:

A v^{i+1} = b

where

A = M + t C + t^2 K b = M v^i - t (K x^i + f_0 + f_plas - f_ext)

In this implementation, the following holds:

M is a diagonal matrix with diagonal elements given by m_mass. C is a diagonal matrix given by massCoef*M, which is Raleigh damping with stiffness coefficient zero.

Also plastic forces, f_plas, is ignored.

Parameters:

	dt	The time step, t, which is about to be taken.
	mass_damping	Coefficient for mass damping in the Raleigh damping equation. The coefficient in C = M + K. In this implementation = 0.

template<typename tetrahedron_iterator , typename real_type >

void OpenTissue::fem::detail::initialize_plastic	(	tetrahedron_iterator	begin,
		tetrahedron_iterator	end,
		real_type const &	c_yield,
		real_type const &	c_creep,
		real_type const &	c_max
	)			[inline]

Initialize plasticity material parameters.

Note: Initialize_stiffness_elements must have been invoked prior to calling this function.

Parameters:

	begin	Iterator to first tetrahedron.
	end	Iterator to one past last tetrahedron.
	c_yield	Plastic yield.
	c_creep	Plastic creep.
	c_max	Plastic max.

template<typename tetrahedron_iterator >

void OpenTissue::fem::detail::initialize_stiffness_elements	(	tetrahedron_iterator	begin,
		tetrahedron_iterator	end
	)			[inline]

NOTE: Material parameters (young modulus and poisson ratio) must have been set prior to invoking this function.

Parameters:

	begin
	end

template<typename fem_mesh , typename real_type >

void OpenTissue::fem::detail::position_update	(	fem_mesh &	mesh,
		real_type const &	dt
	)			[inline]

Position Update.

Note: Velocities must have been computed prior to invoking this function.

Parameters:

	mesh
	dt

template<typename tetrahedron_iterator >

void OpenTissue::fem::detail::reset_orientation	(	tetrahedron_iterator const &	begin,
		tetrahedron_iterator const &	end
	)			[inline]

template<typename tetrahedron_iterator >

void OpenTissue::fem::detail::stiffness_assembly	(	tetrahedron_iterator const &	begin,
		tetrahedron_iterator const &	end
	)			[inline]

Stiffness Matrix Assembly. Observe that prior to invocations of this method the rotations of all tetrahedra must have been setted. Also all stiffness assembly data should have been cleared prior to invocations by using the clear_stiffness_assembly method.

From linear elastostatics we have

K u = f

where u is vertex displacements and f is the resulting nodal forces and K is the stiffness matrix. Using u = x - x0, where x is current position and x0 original position we have

K (x - x0) = f K x - K x0 = f

Introducing f0 = - K x0 as the force offset vectors we have

f = K x + f0

The idea of stiffness warping is to rotate the x-position back into the original coordinate system in which K were original computed, then compute nodal forces in this frame and finally rotate teh nodal forces back to the current frame.

Let R denote the current orientation, then R^{-1} rotates back to the orignal frame which means we have

f = R K ( R^{-1} x - x0 ) f = R K R^{-1} x - R K x0

So with stiffness warping we have to compute

K' = R K R^{-1} f0' = - R K x0

For n nodes the system stiffness matrix K' is a 3n X 3n symmetric and sparse matrix. It would be insane to actual allocated such a matrix instead the matrix is stored inside the nodes of the volume mesh.

Each node stores a row of K' and the correspond entries of the f0' vector.

That is the i'th node stores all non-zero 3-by-3 sub-block of the i'th row of the stiffness matrix and it also stores the i'th 3-dimensional subvector of f0

K' f0' j

• . - - - | . | | | | . | | | i-> |......X....| |x| | . | | | | . | | |
• . - - -

Let M denote the set of all tetrahedral elements then the assembly of the warped stiffness matrix can be written as

K'_ij = sum Re Ke_ij Re^T e M and i, j e

And the assembly of the force offset vector

f0'_i = sum - Re Ke_ij x0_j e M and i, j e

Notice that this can be optimized if node i as a meber of the e'th element then the contribution to the above summation can be written as (recal the indices j,k and m denote the three ofther nodes of e):

-Re Ke_ii x0_i - Re Ke_ij x0_j -Re Ke_im x0_m - Re Ke_ik x0_k -Re ( Ke_ii x0_i + Ke_ij x0_j + Ke_im x0_m + Ke_ik x0_k)

which saves us a few matrix multiplications and we then finally have

f0'_i = sum -Re ( Ke_ii x0_i + Ke_ij x0_j + Ke_im x0_m + Ke_ik x0_k) e M and i e

Assuming that f0'_i is initially cleared to zero for all i. This result in the implementation strategy:

for each element e do for each node i of e do tmp = 0 for each node j of e do tmp += Ke_ij * xo_j next j f0'_i -= Re*tmp next i next e

Also the stiffness matrix can be optimized slightly by exploiting the symmetry property. The symmetry indicates that it is sufficient to only compute the upper triangular and diagonal parts. The lower triangular parts can be found simply be taking the transpose of the upper triangular parts.

So the e'th tetrahedral element contributes with

K'_ij += Re Ke_ij Re^{T} K'_ji += Re Ke_ji Re^{T} = ( Re Ke_ij Re^{T} )^T

because Ke_ji = Ke_ij^T. Assuming that all K'_ij is initially cleared to zero this results in the implementation strategy:

for each element e do for each node i of e do for each node j of e do if i >= j then tmp = Re Ke_ij Re^{T} K'_ij += tmp if j > i then K'_ji += trans(tmp) end if end if next j next i next e

The two implementation strategies can now be combined into one, which is what we have implemented below.

Note that if Re is initially set to the identity matrix, then the stiffness warping reduces to the traditional assembly of stiffness matrix (as it is used in linear elastostatics).

If the i'th node is set to be fixed, this actually corresponds to letting K'_ii = identity and letting K'_ij and K'_ji be zero for all j not equal to i. This is known as a Direchlet boundary condition. However, we do not this during our assembly. Instead we assemble the K' matrix as though there were no fixed nodes. Later on when we use the K' matrix in computations such as K' x, we simply test whether x_j is fixed when it is multilpied by the j'th column of K' i.e. K'_*j. If so we simply do nothing!!! This is computatonally more tracjktable and it also allow us to more easily turn nodes fixed and un-fixed dynamically during animation.

Parameters:

	begin
	end

template<typename real_type , typename tetrahedron_iterator >

void OpenTissue::fem::detail::uniform_density	(	tetrahedron_iterator const &	begin,
		tetrahedron_iterator const &	end,
		real_type const &	density
	)			[inline]

Set Uniform Density.

Parameters:

	begin
	end
	density

template<typename real_type , typename tetrahedron_iterator >

void OpenTissue::fem::detail::uniform_poisson	(	tetrahedron_iterator const &	begin,
		tetrahedron_iterator const &	end,
		real_type const &	poisson
	)			[inline]

Parameters:

	begin
	end
	poisson

template<typename real_type , typename tetrahedron_iterator >

void OpenTissue::fem::detail::uniform_young	(	tetrahedron_iterator const &	begin,
		tetrahedron_iterator const &	end,
		real_type const &	young
	)			[inline]

Parameters:

	begin
	end
	young

template<typename tetrahedron_iterator >

void OpenTissue::fem::detail::update_orientation	(	tetrahedron_iterator const &	begin,
		tetrahedron_iterator const &	end
	)			[inline]