add btFileUtils::toLower to convert a string (char*) to lower case
URDF import demo: add COLLADA .dae file support add FiniteElementMethod demo, extracted from the OpenTissue library (under the zlib license) don't crash if loading an invalid STL file add comparison with Assimp for COLLADA file loading (disabled by default, to avoid library dependency) Gwen: disable some flags that make the build incompatible
This commit is contained in:
erwin coumans
17
Demos3/FiniteElementMethod/OpenTissue/dynamics/fem/fem.h
Normal file
17
Demos3/FiniteElementMethod/OpenTissue/dynamics/fem/fem.h
Normal file
@@ -0,0 +1,17 @@
|
||||
#ifndef OPENTISSUE_DYNAMICS_FEM_FEM_H
|
||||
#define OPENTISSUE_DYNAMICS_FEM_FEM_H
|
||||
//
|
||||
// OpenTissue Template Library
|
||||
// - A generic toolbox for physics-based modeling and simulation.
|
||||
// Copyright (C) 2008 Department of Computer Science, University of Copenhagen.
|
||||
//
|
||||
// OTTL is licensed under zlib: http://opensource.org/licenses/zlib-license.php
|
||||
//
|
||||
#include <OpenTissue/configuration.h>
|
||||
|
||||
//#include <OpenTissue/dynamics/fem/fem_mesh.h>
|
||||
#include <OpenTissue/dynamics/fem/fem_init.h>
|
||||
#include <OpenTissue/dynamics/fem/fem_simulate.h>
|
||||
|
||||
//OPENTISSUE_DYNAMICS_FEM_FEM_H
|
||||
#endif
|
||||
@@ -0,0 +1,221 @@
|
||||
#ifndef OPENTISSUE_DYNAMICS_FEM_FEM_ADD_PLASTICITY_FORCE_H
|
||||
#define OPENTISSUE_DYNAMICS_FEM_FEM_ADD_PLASTICITY_FORCE_H
|
||||
//
|
||||
// OpenTissue Template Library
|
||||
// - A generic toolbox for physics-based modeling and simulation.
|
||||
// Copyright (C) 2008 Department of Computer Science, University of Copenhagen.
|
||||
//
|
||||
// OTTL is licensed under zlib: http://opensource.org/licenses/zlib-license.php
|
||||
//
|
||||
#include <OpenTissue/configuration.h>
|
||||
|
||||
namespace OpenTissue
|
||||
{
|
||||
namespace fem
|
||||
{
|
||||
namespace detail
|
||||
{
|
||||
|
||||
|
||||
template <typename tetrahedron_type, typename real_type >
|
||||
void add_plasticity_force_single1(tetrahedron_type& tet, real_type const & dt)
|
||||
{
|
||||
|
||||
using std::min;
|
||||
using std::sqrt;
|
||||
|
||||
typedef typename tetrahedron_type::vector3_type vector3_type;
|
||||
|
||||
|
||||
assert(tet.m_yield>=0 || !"add_plasticity_force(): yield must be non-negative");
|
||||
assert(tet.m_creep>=0 || !"add_plasticity_force(): creep must be non-negative");
|
||||
assert(tet.m_creep<=(1.0/dt) || !"add_plasticity_force(): creep must be less that reciprocal time-step");
|
||||
assert(tet.m_max>=0 || !"add_plasticity_force(): max must be non-negative");
|
||||
|
||||
//--- Storage for total and elastic strains (plastic strains are stored in the tetrahedra)
|
||||
real_type e_total[6];
|
||||
real_type e_elastic[6];
|
||||
|
||||
for(int i=0;i<6;++i)
|
||||
e_elastic[i] = e_total[i] = 0;
|
||||
|
||||
//--- Compute total strain: e_total = Be (Re^{-1} x - x0)
|
||||
for(unsigned int j=0;j<4;++j)
|
||||
{
|
||||
vector3_type & x_j = tet.node(j)->m_coord;
|
||||
vector3_type & x0_j = tet.node(j)->m_model_coord;
|
||||
|
||||
vector3_type tmp = (trans(tet.m_Re)*x_j) - x0_j;
|
||||
real_type bj = tet.m_B[j](0);
|
||||
real_type cj = tet.m_B[j](1);
|
||||
real_type dj = tet.m_B[j](2);
|
||||
e_total[0] += bj*tmp(0);
|
||||
e_total[1] += cj*tmp(1);
|
||||
e_total[2] += dj*tmp(2);
|
||||
e_total[3] += cj*tmp(0) + bj*tmp(1);
|
||||
e_total[4] += dj*tmp(0) + bj*tmp(2);
|
||||
e_total[5] += dj*tmp(1) + cj*tmp(2);
|
||||
}
|
||||
|
||||
//--- Compute elastic strain
|
||||
for(int i=0;i<6;++i)
|
||||
e_elastic[i] = e_total[i] - tet.m_plastic[i];
|
||||
|
||||
//--- if elastic strain exceeds c_yield then it is added to plastic strain by c_creep
|
||||
real_type norm_elastic = 0;
|
||||
for(int i=0;i<6;++i)
|
||||
norm_elastic += e_elastic[i]*e_elastic[i];
|
||||
norm_elastic = sqrt(norm_elastic);
|
||||
//max_elastic = max(max_elastic,norm_elastic);
|
||||
|
||||
if(norm_elastic > tet.m_yield)
|
||||
{
|
||||
real_type amount = dt*min(tet.m_creep,(1.0/dt)); //--- make sure creep do not exceed 1/dt
|
||||
for(int i=0;i<6;++i)
|
||||
tet.m_plastic[i] += amount*e_elastic[i];
|
||||
}
|
||||
|
||||
//--- if plastic strain exceeds c_max then it is clamped to maximum magnitude
|
||||
real_type norm_plastic = 0;
|
||||
for(int i=0;i<6;++i)
|
||||
norm_plastic += tet.m_plastic[i]*tet.m_plastic[i];
|
||||
norm_plastic = sqrt(norm_plastic);
|
||||
//max_plastic = max(max_plastic,norm_plastic);
|
||||
|
||||
if(norm_plastic > tet.m_max)
|
||||
{
|
||||
real_type scale = tet.m_max/norm_plastic;
|
||||
for(int i=0;i<6;++i)
|
||||
tet.m_plastic[i] *= scale;
|
||||
}
|
||||
//--- Compute plastic forces: f_plastic = Re Pe e_plastic; where Pe = Ve Be^T E
|
||||
for(unsigned int j=0;j<4;++j)
|
||||
{
|
||||
real_type * plastic = tet.m_plastic;
|
||||
real_type bj = tet.m_B[j](0);
|
||||
real_type cj = tet.m_B[j](1);
|
||||
real_type dj = tet.m_B[j](2);
|
||||
real_type E0 = tet.m_D(0);
|
||||
real_type E1 = tet.m_D(1);
|
||||
real_type E2 = tet.m_D(2);
|
||||
vector3_type f;
|
||||
|
||||
//---
|
||||
//---
|
||||
//--- Recall the structure of the B and E matrices
|
||||
//---
|
||||
//--- | bj 0 0 | | E0 E1 E1 |
|
||||
//--- B_j = | 0 cj 0 | | E1 E0 E0 |
|
||||
//--- | 0 0 dj | E = | E1 E1 E0 |
|
||||
//--- | cj bj 0 | | E2 |
|
||||
//--- | dj 0 bj | | E2 |
|
||||
//--- | 0 dj cj | | E2 |
|
||||
//---
|
||||
//--- This implyies that the product B_j^T E is
|
||||
//---
|
||||
//--- | bj E0 bj E1 bj E1 cj E2 dj E2 0 |
|
||||
//--- | cj E1 cj E0 cj E1 bj E2 0 dj E2 |
|
||||
//--- | dj E1 dj E1 dj E0 0 bj E2 cj E2 |
|
||||
//---
|
||||
//--- Notice that eventhough this is a 3X6 matrix with 18-3=15 nonzero elements
|
||||
//--- it actually only contains 9 different value.
|
||||
//---
|
||||
//--- In fact these values could be precomputed and stored in each tetrahedron.
|
||||
//--- Furthermore they could be pre-multiplied by the volume of the tetrahedron
|
||||
//--- to yield the matrix,
|
||||
//---
|
||||
//--- P_j = Ve B_j^T E
|
||||
//---
|
||||
//--- The plastic force that should be subtracted from node j is then computed
|
||||
//--- in each iteration as
|
||||
//---
|
||||
//--- f_j = Re P_j e_plastic
|
||||
//---
|
||||
//---
|
||||
real_type bjE0 = bj*E0;
|
||||
real_type bjE1 = bj*E1;
|
||||
real_type bjE2 = bj*E2;
|
||||
real_type cjE0 = cj*E0;
|
||||
real_type cjE1 = cj*E1;
|
||||
real_type cjE2 = cj*E2;
|
||||
real_type djE0 = dj*E0;
|
||||
real_type djE1 = dj*E1;
|
||||
real_type djE2 = dj*E2;
|
||||
|
||||
f(0) = bjE0*plastic[0] + bjE1*plastic[1] + bjE1*plastic[2] + cjE2*plastic[3] + djE2*plastic[4];
|
||||
f(1) = cjE1*plastic[0] + cjE0*plastic[1] + cjE1*plastic[2] + bjE2*plastic[3] + + djE2*plastic[5];
|
||||
f(2) = djE1*plastic[0] + djE1*plastic[1] + djE0*plastic[2] + bjE2*plastic[4] + cjE2*plastic[5];
|
||||
|
||||
f *= tet.m_V;
|
||||
tet.node(j)->m_f_external += tet.m_Re*f;
|
||||
}
|
||||
}
|
||||
/**
|
||||
* Add plasticity forces.
|
||||
*
|
||||
* @param begin Iterator to first tetrahedron.
|
||||
* @param end Iterator to one past last tetrahedron.
|
||||
* @param dt Simulation time step
|
||||
*
|
||||
*/
|
||||
template < typename tetrahedron_iterator,typename real_type >
|
||||
inline void add_plasticity_force2(
|
||||
tetrahedron_iterator begin
|
||||
, tetrahedron_iterator end
|
||||
, real_type const & dt
|
||||
)
|
||||
{
|
||||
//typedef typename tetrahedron_iterator::value_type::matrix3x3_type matrix3x3_type;
|
||||
|
||||
//real_type max_elastic = 0;
|
||||
//real_type max_plastic = 0;
|
||||
|
||||
assert(dt>0 || !"add_plasticity_force(): time-step must be positive");
|
||||
//---
|
||||
//--- The total strain is given by
|
||||
//---
|
||||
//--- e_total = Be u
|
||||
//---
|
||||
//--- Where u is the nodal displacement, ie u = x - x0. Applying
|
||||
//--- the idea of stiffness warping we have
|
||||
//---
|
||||
//--- e_total = Be ( Re^{-1} x - x0)
|
||||
//---
|
||||
//--- where R_e is the rotational deformation of the e'th
|
||||
//--- tetrahedral element.
|
||||
//---
|
||||
//--- The elastic strain is given as the total strain minus the plastic strain
|
||||
//---
|
||||
//--- e_elastic = e_total - e_plastic
|
||||
//---
|
||||
//--- e_plastic is initial initialized to zero. Plastic strain causes plastic forces in the material
|
||||
//---
|
||||
//--- f_plastic = Re Ke u_plastic
|
||||
//---
|
||||
//--- Using e_plastic = Be u_plastic
|
||||
//---
|
||||
//--- f_plastic = Re Ke Be^{-1} e_plastic
|
||||
//--- f_plastic = Re (Ve Be^T E Be ) Be^{-1} e_plastic
|
||||
//--- f_plastic = Re (Ve Be^T E) e_plastic
|
||||
//---
|
||||
//--- Introducing the plasticity matrix Pe = (Ve Be^T E) we have
|
||||
//---
|
||||
//--- f_plastic = Re Pe e_plastic
|
||||
//---
|
||||
|
||||
for (tetrahedron_iterator T = begin; T != end; ++T)
|
||||
{
|
||||
add_plasticity_force_single(T,dt);
|
||||
|
||||
}
|
||||
//std::cout << "max elastic = " << max_elastic << std::endl;
|
||||
//std::cout << "max plastic = " << max_plastic << std::endl;
|
||||
|
||||
}
|
||||
|
||||
} // namespace detail
|
||||
} // namespace fem
|
||||
} // namespace OpenTissue
|
||||
|
||||
//OPENTISSUE_DYNAMICS_FEM_FEM_ADD_PLASTICITY_FORCE_H
|
||||
#endif
|
||||
@@ -0,0 +1,37 @@
|
||||
#ifndef OPENTISSUE_DYNAMICS_FEM_FEM_CLEAR_STIFFNESS_ASSEMBLY_H
|
||||
#define OPENTISSUE_DYNAMICS_FEM_FEM_CLEAR_STIFFNESS_ASSEMBLY_H
|
||||
//
|
||||
// OpenTissue Template Library
|
||||
// - A generic toolbox for physics-based modeling and simulation.
|
||||
// Copyright (C) 2008 Department of Computer Science, University of Copenhagen.
|
||||
//
|
||||
// OTTL is licensed under zlib: http://opensource.org/licenses/zlib-license.php
|
||||
//
|
||||
#include <OpenTissue/configuration.h>
|
||||
|
||||
namespace OpenTissue
|
||||
{
|
||||
namespace fem
|
||||
{
|
||||
namespace detail
|
||||
{
|
||||
/**
|
||||
*
|
||||
* @param begin
|
||||
* @param end
|
||||
*/
|
||||
template < typename node_type >
|
||||
inline void clear_stiffness_assembly_single(node_type* node)
|
||||
{
|
||||
typedef typename node_type::matrix_iterator matrix_iterator;
|
||||
node->m_f0.clear();
|
||||
for (matrix_iterator Kij = node->Kbegin() ; Kij != node->Kend(); ++Kij )
|
||||
Kij->second.clear();
|
||||
}
|
||||
|
||||
} // namespace detail
|
||||
} // namespace fem
|
||||
} // namespace OpenTissue
|
||||
|
||||
//OPENTISSUE_DYNAMICS_FEM_FEM_CLEAR_STIFFNESS_ASSEMBLY_H
|
||||
#endif
|
||||
@@ -0,0 +1,403 @@
|
||||
#ifndef OPENTISSUE_DYNAMICS_FEM_FEM_COMPUTE_B_H
|
||||
#define OPENTISSUE_DYNAMICS_FEM_FEM_COMPUTE_B_H
|
||||
//
|
||||
// OpenTissue Template Library
|
||||
// - A generic toolbox for physics-based modeling and simulation.
|
||||
// Copyright (C) 2008 Department of Computer Science, University of Copenhagen.
|
||||
//
|
||||
// OTTL is licensed under zlib: http://opensource.org/licenses/zlib-license.php
|
||||
//
|
||||
#include <OpenTissue/configuration.h>
|
||||
|
||||
namespace OpenTissue
|
||||
{
|
||||
namespace fem
|
||||
{
|
||||
namespace detail
|
||||
{
|
||||
|
||||
/**
|
||||
* Calculate B-matrices (derivatives of shape functions N, i.e B = SN).
|
||||
*
|
||||
* @param e10 Given the four corner points p0,p1,p2, and p3. This is p1-p0.
|
||||
* @param e20 Given the four corner points p0,p1,p2, and p3. This is p2-p0.
|
||||
* @param e30 Given the four corner points p0,p1,p2, and p3. This is p3-p0.
|
||||
* @param volume The volume of the tetrahedron.
|
||||
* @param B Upon return this array of 4 B-vectors contains the computed values.
|
||||
*/
|
||||
template<typename real_type, typename vector3_type>
|
||||
inline void compute_B(
|
||||
vector3_type const& e10,
|
||||
vector3_type const& e20,
|
||||
vector3_type const& e30,
|
||||
real_type volume,
|
||||
vector3_type * B
|
||||
)
|
||||
{
|
||||
//
|
||||
// In summary we want to calculate
|
||||
//
|
||||
// B = S N
|
||||
//
|
||||
// where
|
||||
//
|
||||
// B = [ B0 B1 B2 B3 ], N = [ N0 N1 N2 N3 ]
|
||||
//
|
||||
// Here S is the operational matrix given by
|
||||
// - -
|
||||
// | d/dx 0 0 |
|
||||
// | 0 d/dy 0 |
|
||||
// | 0 0 d/dz |
|
||||
// S = | d/dy d/dx 0 |
|
||||
// | d/dz 0 d/dx |
|
||||
// | 0 d/dz d/dy |
|
||||
// - -
|
||||
// N_i's are called the shape functions and they are used to interpolate values. For
|
||||
// 3D linear elastostatics these can be written as
|
||||
//
|
||||
// - -
|
||||
// | w_i 0 0 |
|
||||
// N_i = | 0 w_i 0 |
|
||||
// | 0 0 w_i |
|
||||
// - -
|
||||
//
|
||||
// The w_i's are functions in the coordinates x,y, and z and will be derived
|
||||
// shortly below of here.
|
||||
//
|
||||
// The matrix B_i is given by:
|
||||
//
|
||||
// | b_i 0 0 |
|
||||
// | 0 c_i 0 |
|
||||
// B_i = | 0 0 d_i |
|
||||
// | c_i b_i 0 |
|
||||
// | d_i 0 b_i |
|
||||
// | 0 d_i c_i |
|
||||
//
|
||||
// The entries are stored as
|
||||
//
|
||||
// B[i] = [ b_i c_i d_i ]
|
||||
//
|
||||
// In the following we go into details of deriving formulas for the interpolating
|
||||
// functions N_i. From these formulas it follows how to compute the derivatives.
|
||||
// Finally we show a clever way of implementing the computaitons.
|
||||
//
|
||||
// Using a linear polynomial for interpolating a u(x,y,z)-value at a given
|
||||
// position x,y,z inside the tetrahedron
|
||||
//
|
||||
// u(x,y,z) = a_0 + a_1 x + a_2 y + a_3 z = P^T A
|
||||
//
|
||||
// where the a's are the polynomial coefficients (to be determed) and
|
||||
//
|
||||
// | 1 | | a0 |
|
||||
// P = | x | , and A = | a1 |
|
||||
// | y | | a2 |
|
||||
// | z | | a3 |
|
||||
//
|
||||
// Say we know the u-values, u0(x0,y0,z0),...,u3(x3,y2,z3), at the four tetrahedron corner
|
||||
// points, P0^T [x0 y0 z0],...,P3^T=[x3 y3 z3], then we can set up the linear system
|
||||
//
|
||||
// | u0 | | 1 x0 y0 z0 | | a0 |
|
||||
// | u1 | = | 1 x1 y1 z1 | | a1 | = C A
|
||||
// | u2 | | 1 x2 y2 z2 | | a2 |
|
||||
// | u3 | | 1 x3 y3 z3 | | a3 |
|
||||
//
|
||||
// The C matrix is invertible (as long as the four points are in general postion, meaning
|
||||
// that no point can be written as a linear combination of any of the three other points),
|
||||
// so we can solve the polynomial coefficients by
|
||||
//
|
||||
// | u0 |
|
||||
// A = C^-1 | u1 |
|
||||
// | u2 |
|
||||
// | u3 |
|
||||
// Substitution into the polynomial interpolation formula yields
|
||||
//
|
||||
// | u0 |
|
||||
// u(x,y,z) = P^T A = P^T C^{-1} | u1 | (*1)
|
||||
// | u2 |
|
||||
// | u3 |
|
||||
//
|
||||
// Denoting the i'th column of, P^T C^{-1}, by N_i we have u(x,y,z) = sum_i N_i u_i
|
||||
//
|
||||
//
|
||||
// Let us now try to look at the bary-centric coordinates, w0,...,w1 of (x,y,z), these
|
||||
// can also be used for interpolation
|
||||
//
|
||||
// | u0 | |u0|
|
||||
// u(z,y,z) = w0 u0 + w1 u1 + w2 u2 + w3 u3 = [w0 w1 w2 w3] | u1 | = W^T |u1|
|
||||
// | u2 | |u2|
|
||||
// | u3 | |u3|
|
||||
//
|
||||
// From the coordinates of the four points and the condition 1 = w0 + w1 + w2 + w3 we can set
|
||||
// up the linear system
|
||||
//
|
||||
// | 1 | | 1 1 1 1 | | w0 |
|
||||
// | x | = | x0 x1 x2 x3 | | w1 |
|
||||
// | y | | y0 y1 y2 y3 | | w2 |
|
||||
// | z | | z0 z1 z2 z3 | | w3 |
|
||||
//
|
||||
//
|
||||
// P = Q W
|
||||
//
|
||||
// Since the four points are in general postion, Q is invertible and we can solve for W
|
||||
//
|
||||
// W = Q^{-1} P
|
||||
//
|
||||
// Insertion into the barycentric interpolation formula yields
|
||||
//
|
||||
// |u0| |u0| |u0|
|
||||
// u(x,y,z) = W^T |u1| = ( Q^{-1} P )^T |u1| = P^T Q^{-T} |u1| (*2)
|
||||
// |u2| |u2| |u2|
|
||||
// |u3| |u3| |u3|
|
||||
//
|
||||
// Comparision with the polynomial interpolation derivation we see that C = Q^T, furthermore
|
||||
// we notice that w_i = N_i. So (not very surprinsingly) bary-centric interpolation is really
|
||||
// just linear polynomial interpolation.
|
||||
//
|
||||
// Notice that for the volume, V, of the tetrahedron we have the relation: 6 V = det( C ) = det ( Q^T )
|
||||
//
|
||||
// When computing the stiffness matrix we are interested in derivatives of the N_i's with
|
||||
// respect to the x,y and z coordinates.
|
||||
//
|
||||
//
|
||||
// From (*1) and (*2) we see (recall we use zero-indexing)
|
||||
//
|
||||
// d N_i
|
||||
// ----- = C^{-1}_{1,i} = Q^{-T}_{1,i} (*3a)
|
||||
// d x
|
||||
//
|
||||
// d N_i
|
||||
// ----- = C^{-1}_{2,i} = Q^{-T}_{2,i} (*3b)
|
||||
// d y
|
||||
//
|
||||
// d N_i
|
||||
// ----- = C^{-1}_{3,i} = Q^{-T}_{3,i} (*3c)
|
||||
// d z
|
||||
//
|
||||
// Instead of actually computing the inverse of the 4x4 C or Q matrices a computational
|
||||
// more efficient solution can be derived, which only requires us to solve for a 3x3 system.
|
||||
//
|
||||
// By Cramers rule we have
|
||||
//
|
||||
// (-1)^{i+j} det( C_ji)
|
||||
// C^{-1}_{i,j} = ----------------------
|
||||
// det(C)
|
||||
//
|
||||
// where det(C_ji) is the determinant of C with the j'th row and i'th column removed.
|
||||
//
|
||||
// Now defined
|
||||
//
|
||||
// e10 = P1 - P0
|
||||
// e20 = P2 - P0
|
||||
// e30 = P3 - P0
|
||||
//
|
||||
// and the matrix E
|
||||
//
|
||||
// | e10^T | | (x1-x0) (y1-y0) (z1-z0) |
|
||||
// E = | e20^T | = | (x2-x0) (y2-y0) (z2-z0) |
|
||||
// | e30^T | | (x3-x0) (y3-y0) (z3-z0) |
|
||||
//
|
||||
// Then
|
||||
//
|
||||
// det(C) = det(E) and C^{-1}_{i+1,j+1} = E^{-1}_{i,j} (*4)
|
||||
//
|
||||
// This can shown by straightforward computation, immediately is is verified that
|
||||
//
|
||||
// {-1}^((i+1)+(j+1)) = {-1}^(i+j)
|
||||
//
|
||||
// So we have to show
|
||||
//
|
||||
// det( C_(i+1,j+1)) det( E_(ij))
|
||||
// --------------------- = --------------
|
||||
// det(C) det(E)
|
||||
//
|
||||
// assume det(C)= det(E) (left for reader as exercise:-) then inorder to prove
|
||||
// second half of (*4) we have to show
|
||||
//
|
||||
// det( C_(i+1,j+1)) = det( E_(ij))
|
||||
//
|
||||
// A total of 9 cases exist, for here we will show a single case and leave
|
||||
// the remaining cases for the reader, use i=1 and j= 0, implying
|
||||
//
|
||||
// det( C_(12)) = det( E_(01))
|
||||
//
|
||||
// So
|
||||
// - -
|
||||
// | 1 x0 z0 | - -
|
||||
// det | 1 x2 z2 | = det | (x2-x0) (z2-z0)|
|
||||
// | 1 x3 z3 | | (x3-x0) (z3-z0)|
|
||||
// - - - -
|
||||
//
|
||||
// Which is trivially true. So we have
|
||||
//
|
||||
// | . ... | | . ... |
|
||||
// C^{-1} = | . E^{-1} | ; Q^{-T} = | . E^{-T} |
|
||||
//
|
||||
// As can be seen from (*3) we do not have all the values needed since
|
||||
// the first columns of C^{-1} and Q^{-T} are missing. To remedy this problem we can use
|
||||
// the condition of the bary-centric coordinates
|
||||
//
|
||||
// w0 = 1 - w1 - w2 - w3
|
||||
//
|
||||
// Taking the derivate wrt. x,y, and z yields
|
||||
//
|
||||
// d N_0
|
||||
// ----- = 1 - E^{-1}_01 - E^{-1}_02 - E^{-1}_03 (*5a)
|
||||
// d x
|
||||
//
|
||||
// d N_0
|
||||
// ----- = 1 - E^{-1}_11 - E^{-1}_12 - E^{-1}_13 (*5b)
|
||||
// d y
|
||||
//
|
||||
// d N_0
|
||||
// ----- = 1 - E^{-1}_21 - E^{-1}_22 - E^{-1}_23 (*5c)
|
||||
// d z
|
||||
//
|
||||
//
|
||||
// and for i>0 we have
|
||||
//
|
||||
// d N_i
|
||||
// ----- = E^{-1}_0i (*6a)
|
||||
// d x
|
||||
//
|
||||
// d N_i
|
||||
// ----- = E^{-1}_1i (*6b)
|
||||
// d y
|
||||
//
|
||||
//
|
||||
// d N_i
|
||||
// ----- = E^{-1}_2i (*6c)
|
||||
// d z
|
||||
//
|
||||
// The notation b_i = d N_i / dx, c_i = d N_i / dy, and d_i = d N_i / dz is used. Furthermore
|
||||
// all the derivatives are returned as four B-vectors, where
|
||||
//
|
||||
// B[i] = [ b_i c_i d_i]
|
||||
//
|
||||
//
|
||||
//
|
||||
// The above may seem as a mathematical trick, so let us try to derive our
|
||||
// equations a litlle differently, but before doing so we must first develop
|
||||
// some equations for the volume of a tetrahedron.
|
||||
//
|
||||
// | e10^T |
|
||||
// 6 V = det(E) = | e20^T | = e10 \cdot (e20 \times e30)
|
||||
// | e30^T |
|
||||
//
|
||||
// where e10 = p1-p0, e20 = p2-p0 and so on.
|
||||
// In general given the right-hand order i,j,k, and m of the nodes we
|
||||
// write vol(i,j,k,m) = eji \cdot (eki \times emi) / 6
|
||||
//
|
||||
// Barycentric coordinates are infact the volume weighted weights coresponding to
|
||||
// the four tetrahedra lying inside the tetrahedron and having apex at the
|
||||
// point p=[x,y,z]^T and bases equal to the 4 triangular faces of the enclosing
|
||||
// tetrahedron. To realize this we will examine bary-centric coordinates
|
||||
// in the 2D case, that is the case of the Triangle. In 2D area corresponds
|
||||
// to the volume, so given a trianle and a point p inside it
|
||||
//
|
||||
// + 2
|
||||
// /|
|
||||
// / |
|
||||
// / |
|
||||
// / |
|
||||
// / + p|
|
||||
// / |
|
||||
// / |
|
||||
// 0 +-------+ 1
|
||||
//
|
||||
// Let the area weighted weights be defined as
|
||||
//
|
||||
// w_0 = area(1,2,p) / area(0,1,2)
|
||||
// w_1 = area(0,P,2) / area(0,1,2) = area(2,0,P) / area(0,1,2)
|
||||
// w_2 = area(0,1,p) / area(0,1,2)
|
||||
//
|
||||
// It is intuitive to see that if p moves towards the i'th corner point then
|
||||
// the nominator of w_i goes towards area(0,1,2). This means that w_i -> 1
|
||||
// while w_j ->0 for j\neq i.
|
||||
//
|
||||
// Having introduced the barycentric coordinates as the area weighted weights
|
||||
// it is quite intuitive to see that we also must have
|
||||
//
|
||||
// 1 = \sum_i w_i
|
||||
//
|
||||
// In 3D the barycentric coordinates are the volume weighted weights. That is
|
||||
// given the right handed order 0,1,2 and 3 of the corner points, we have by
|
||||
// analogy to the 2D case
|
||||
//
|
||||
// vol(0,1,2,p)
|
||||
// w_3 = ------------ = (p-p0) \codt ( e10 \times \e20 ) / 6V
|
||||
// V
|
||||
//
|
||||
// vol(0,1,3,p)
|
||||
// w_2 = ------------ = (p-p0) \codt ( e10 \times \e30 ) / 6V
|
||||
// V
|
||||
//
|
||||
// vol(0,2,3,p)
|
||||
// w_1 = ------------ = (p-p0) \codt ( e20 \times \e30 ) / 6V
|
||||
// V
|
||||
//
|
||||
// vol(1,3,2,p)
|
||||
// w_0 = ------------ = (p-p0) \codt ( e31 \times \e21 ) / 6V
|
||||
// V
|
||||
//
|
||||
// But since we allready know w_3, w_2 and w_1 it is faster to compute w_0 as
|
||||
//
|
||||
// w_0 = 1 - w_1 - w_2 - w_3
|
||||
//
|
||||
// Recal from previous that we seek the derivatives of w_i's wrt. the coordinaes when
|
||||
// we compute the B functions. Let us write the [x y z] coordinates as [x_0 x_1 x_2] then
|
||||
//
|
||||
// d w_3
|
||||
// ------- = ( e10 \times \e20 )_{x_j} / 6V
|
||||
// d x_j
|
||||
//
|
||||
// That is the j'th coordinate of the cross product divided by 6 times the volume. Similar for
|
||||
//
|
||||
// d w_2
|
||||
// ------- = ( e10 \times \e30 )_{x_j} / 6V
|
||||
// d x_j
|
||||
//
|
||||
// d w_1
|
||||
// ------- = ( e20 \times \e30 )_{x_j} / 6V
|
||||
// d x_j
|
||||
//
|
||||
// and finally
|
||||
//
|
||||
// d w_0 d w_i
|
||||
// ------- = 1 - sum_i>0 -------
|
||||
// d x_j d x_j
|
||||
//
|
||||
// Our formulas may seem differnt from those in (*6), however
|
||||
//
|
||||
// d w_k
|
||||
// ------- = ( eji \times \emi )_{x_j} / 6V = E_{-1}_{k,j}
|
||||
// d x_j
|
||||
//
|
||||
// That is the adjoint of E_(j,k) is given by ( eji \times \emi )_{x_j} / 6.
|
||||
//
|
||||
//
|
||||
//
|
||||
|
||||
real_type div6V = 1/(6.0*volume);
|
||||
|
||||
B[1](0) = (e20(2) * e30(1) - e20(1) * e30(2)) * div6V; // b0 = -det(E_11), (where E_ij is the submatrix of E with row i and column j removed)
|
||||
B[2](0) = (e10(1) * e30(2) - e10(2) * e30(1)) * div6V; // b1 = det(E_12)
|
||||
B[3](0) = (e10(2) * e20(1) - e10(1) * e20(2)) * div6V; // b2 = -det(E_13)
|
||||
B[0](0) = -B[1](0) - B[2](0) - B[3](0); // b3 = -b0 - b1 - b2
|
||||
|
||||
B[1](1) = (e20(0) * e30(2) - e20(2) * e30(0)) * div6V; // c0 = det(E_21)
|
||||
B[2](1) = (e10(2) * e30(0) - e10(0) * e30(2)) * div6V; // c1 = -det(E_22)
|
||||
B[3](1) = (e10(0) * e20(2) - e10(2) * e20(0)) * div6V; // c2 = det(E_23)
|
||||
B[0](1) = -B[1](1) - B[2](1) - B[3](1); // c3 = -c0 - c1 - c2
|
||||
|
||||
B[1](2) = (e20(1) * e30(0) - e20(0) * e30(1)) * div6V; // d0 = -det(E_31)
|
||||
B[2](2) = (e10(0) * e30(1) - e10(1) * e30(0)) * div6V; // d1 = det(E_32)
|
||||
B[3](2) = (e10(1) * e20(0) - e10(0) * e20(1)) * div6V; // d2 = -det(E_33)
|
||||
B[0](2) = -B[1](2) - B[2](2) - B[3](2); // d3 = -d0 - d1 - d2
|
||||
}
|
||||
|
||||
} // namespace detail
|
||||
} // namespace fem
|
||||
} // namespace OpenTissue
|
||||
|
||||
//OPENTISSUE_DYNAMICS_FEM_FEM_COMPUTE_B_H
|
||||
#endif
|
||||
@@ -0,0 +1,121 @@
|
||||
#ifndef OPENTISSUE_DYNAMICS_FEM_FEM_COMPUTE_ISOTROPIC_ELASTICITY_H
|
||||
#define OPENTISSUE_DYNAMICS_FEM_FEM_COMPUTE_ISOTROPIC_ELASTICITY_H
|
||||
//
|
||||
// OpenTissue Template Library
|
||||
// - A generic toolbox for physics-based modeling and simulation.
|
||||
// Copyright (C) 2008 Department of Computer Science, University of Copenhagen.
|
||||
//
|
||||
// OTTL is licensed under zlib: http://opensource.org/licenses/zlib-license.php
|
||||
//
|
||||
#include <OpenTissue/configuration.h>
|
||||
|
||||
namespace OpenTissue
|
||||
{
|
||||
namespace fem
|
||||
{
|
||||
namespace detail
|
||||
{
|
||||
/**
|
||||
* Compute Isotropic Elasticity Matrix.
|
||||
*
|
||||
* @param young Young Modulus.
|
||||
* @param poisson Poisson ratio.
|
||||
* @param D Upon return holds the elasticity matrix in vector form D = [D0,D1,D2]
|
||||
*/
|
||||
template<typename real_type, typename vector3_type>
|
||||
inline void compute_isotropic_elasticity_vector(
|
||||
real_type const & young
|
||||
, real_type const & poisson
|
||||
, vector3_type & D
|
||||
)
|
||||
{
|
||||
assert(young>0 || !"compute_isotropic_elasticity_vector(): Young modulus must be positive");
|
||||
assert(poisson>0 || !"compute_isotropic_elasticity_vector(): Poisson ratio must be positive");
|
||||
assert(poisson<0.5 || !"compute_isotropic_elasticity_vector(): Poisson ratio must be less than a half");
|
||||
|
||||
// The isotropic Elasticity matrix D is given by
|
||||
// - - - -
|
||||
// | 1-nu nu nu 0 0 0 | | D0 D1 D1 0 0 0 |
|
||||
// young | nu 1-nu nu 0 0 0 | | D1 D0 D1 0 0 0 |
|
||||
// ------------- | nu nu 1-nu 0 0 0 | | D1 D1 D0 0 0 0 |
|
||||
// (1+nu)(1-2nu) | 0 0 0 (1-2nu)/2 0 0 | = | 0 0 0 D2 0 0 |
|
||||
// | 0 0 0 0 (1-2nu)/2 0 | | 0 0 0 0 D2 0 |
|
||||
// | 0 0 0 0 0 (1-2nu)/2 | | 0 0 0 0 0 D2 |
|
||||
// - - - -
|
||||
|
||||
real_type poisson2 = 2.0*poisson;
|
||||
real_type scale = young / ((1.0 + poisson) * (1.0 - poisson2));
|
||||
D(0) = (1.0 - poisson) * scale;
|
||||
D(1) = poisson * scale;
|
||||
D(2) = young / (2.0 + poisson2);
|
||||
}
|
||||
|
||||
/**
|
||||
* Compute Elasticy Matrix.
|
||||
*
|
||||
* @param young Youngs modulus (stiffness)
|
||||
* @param poisson Poissons ratio (compressability)
|
||||
* @param D Upon return holds the isotropoc linear elasticity matrix.
|
||||
*/
|
||||
template<typename real_type,typename matrix_type>
|
||||
inline void compute_isotropic_elasticity_matix(
|
||||
real_type const & young
|
||||
, real_type const & poisson
|
||||
, matrix_type & D
|
||||
)
|
||||
{
|
||||
assert(young>0 || !"compute_isotropic_elasticity_matix(): Young modulus must be positive");
|
||||
assert(poisson>0 || !"compute_isotropic_elasticity_matix(): Poisson ratio must be positive");
|
||||
assert(poisson<0.5 || !"compute_isotropic_elasticity_matix(): Poisson ratio must be less than a half");
|
||||
assert(D.size1()==6 || !"compute_isotropic_elasticity_matix(): D-matrix did not have 6 rows");
|
||||
assert(D.size2()==6 || !"compute_isotropic_elasticity_matix(): D-matrix did not have 6 columns");
|
||||
|
||||
typedef typename matrix_type::value_type value_type;
|
||||
|
||||
value_type lambda = (poisson*young)/(( 1+poisson )*( 1- (2*poisson) )); //--- Lame modulus
|
||||
value_type mu = young/(2*(1+poisson)); //--- Shear Modulus
|
||||
value_type tmp = lambda+(2*mu);
|
||||
|
||||
D(0,0) = tmp; D(0,1) = lambda; D(0,2) = lambda; D(0,3) = 0; D(0,4) = 0; D(0,5) = 0;
|
||||
D(1,0) = lambda; D(1,1) = tmp; D(1,2) = lambda; D(1,3) = 0; D(1,4) = 0; D(1,5) = 0;
|
||||
D(2,0) = lambda; D(2,1) = lambda; D(2,2) = tmp; D(2,3) = 0; D(2,4) = 0; D(2,5) = 0;
|
||||
D(3,0) = 0; D(3,1) = 0; D(3,2) = 0; D(3,3) = mu; D(3,4) = 0; D(3,4) = 0;
|
||||
D(4,0) = 0; D(4,1) = 0; D(4,2) = 0; D(4,3) = 0; D(4,4) = mu; D(4,5) = 0;
|
||||
D(5,0) = 0; D(5,1) = 0; D(5,2) = 0; D(5,3) = 0; D(5,4) = 0; D(5,5) = mu;
|
||||
}
|
||||
|
||||
/**
|
||||
* Alternative representation of isotropic elasticity matrix.
|
||||
*
|
||||
*
|
||||
* @param young Youngs modulus (stiffness)
|
||||
* @param poisson Poissons ratio (compressability)
|
||||
* @param D00
|
||||
* @param D01
|
||||
* @param D33
|
||||
*
|
||||
*/
|
||||
template<typename real_type>
|
||||
inline void compute_isotropic_elasticity_vector(
|
||||
real_type const & young
|
||||
, real_type const & poisson
|
||||
, real_type & D00
|
||||
, real_type & D01
|
||||
, real_type & D33
|
||||
)
|
||||
{
|
||||
assert(young>0 || !"compute_isotropic_elasticity_vector(): Young modulus must be positive");
|
||||
assert(poisson>0 || !"compute_isotropic_elasticity_vector(): Poisson ratio must be positive");
|
||||
assert(poisson<0.5 || !"compute_isotropic_elasticity_vector(): Poisson ratio must be less than a half");
|
||||
|
||||
D01 = (poisson*young)/(( 1 + poisson )*( 1- (2*poisson) )); //--- Lame modulus
|
||||
D33 = young/(2*(1+poisson)); //--- Shear Modulus
|
||||
D00 = D01+(2*D33);
|
||||
}
|
||||
|
||||
} // namespace detail
|
||||
} // namespace fem
|
||||
} // namespace OpenTissue
|
||||
|
||||
//OPENTISSUE_DYNAMICS_FEM_FEM_COMPUTE_ISOTROPIC_ELASTICITY_H
|
||||
#endif
|
||||
@@ -0,0 +1,299 @@
|
||||
#ifndef OPENTISSUE_DYNAMICS_FEM_FEM_COMPUTE_KE_H
|
||||
#define OPENTISSUE_DYNAMICS_FEM_FEM_COMPUTE_KE_H
|
||||
//
|
||||
// OpenTissue Template Library
|
||||
// - A generic toolbox for physics-based modeling and simulation.
|
||||
// Copyright (C) 2008 Department of Computer Science, University of Copenhagen.
|
||||
//
|
||||
// OTTL is licensed under zlib: http://opensource.org/licenses/zlib-license.php
|
||||
//
|
||||
#include <OpenTissue/configuration.h>
|
||||
|
||||
namespace OpenTissue
|
||||
{
|
||||
namespace fem
|
||||
{
|
||||
namespace detail
|
||||
{
|
||||
/**
|
||||
* Compute Stiffness matrix of tetrahedral element.
|
||||
*
|
||||
* @param p0
|
||||
* @param p1
|
||||
* @param p2
|
||||
* @param p3
|
||||
* @param E Youngs modulus.
|
||||
* @param nu Poisson ratio.
|
||||
* @param Ke Upon return contains the computed stiffness values.
|
||||
*/
|
||||
template<typename real_type, typename vector3_type, typename matrix3x3_type>
|
||||
inline 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]
|
||||
)
|
||||
{
|
||||
using std::fabs;
|
||||
|
||||
real_type d = p0(0);
|
||||
real_type d4 = p0(1);
|
||||
real_type d8 = p0(2);
|
||||
real_type d1 = p1(0) - d;
|
||||
real_type d5 = p1(1) - d4;
|
||||
real_type d9 = p1(2) - d8;
|
||||
real_type d2 = p2(0) - d;
|
||||
real_type d6 = p2(1) - d4;
|
||||
real_type d10 = p2(2) - d8;
|
||||
real_type d3 = p3(0) - d;
|
||||
real_type d7 = p3(1) - d4;
|
||||
real_type d11 = p3(2) - d8;
|
||||
real_type d12 = (d1 * d6 * d11 + d2 * d7 * d9 + d3 * d5 * d10) - d1 * d7 * d10 - d2 * d5 * d11 - d3 * d6 * d9;
|
||||
real_type d13 = 1.0 / d12;
|
||||
real_type d14 = fabs(d12 / 6.0);
|
||||
vector3_type B[4];
|
||||
B[0].clear();
|
||||
B[1].clear();
|
||||
B[2].clear();
|
||||
B[3].clear();
|
||||
B[1](0) = (d10 * d7 - d6 * d11) * d13;
|
||||
B[2](0) = (d5 * d11 - d9 * d7) * d13;
|
||||
B[3](0) = (d9 * d6 - d5 * d10) * d13;
|
||||
B[0](0) = -B[1](0) - B[2](0) - B[3](0);
|
||||
B[1](1) = (d2 * d11 - d10 * d3) * d13;
|
||||
B[2](1) = (d9 * d3 - d1 * d11) * d13;
|
||||
B[3](1) = (d1 * d10 - d9 * d2) * d13;
|
||||
B[0](1) = -B[1](1) - B[2](1) - B[3](1);
|
||||
B[1](2) = (d6 * d3 - d2 * d7) * d13;
|
||||
B[2](2) = (d1 * d7 - d5 * d3) * d13;
|
||||
B[3](2) = (d5 * d2 - d1 * d6) * d13;
|
||||
B[0](2) = -B[1](2) - B[2](2) - B[3](2);
|
||||
real_type d15 = E / (1.0 + nu) / (1.0 - 2 * nu);
|
||||
real_type d16 = (1.0 - nu) * d15;
|
||||
real_type d17 = nu * d15;
|
||||
real_type d18 = E / 2 / (1.0 + nu);
|
||||
for(int i = 0; i < 4; ++i)
|
||||
{
|
||||
for(int j = 0; j < 4; ++j)
|
||||
{
|
||||
real_type d19 = B[i](0);
|
||||
real_type d20 = B[i](1);
|
||||
real_type d21 = B[i](2);
|
||||
real_type d22 = B[j](0);
|
||||
real_type d23 = B[j](1);
|
||||
real_type d24 = B[j](2);
|
||||
Ke[i][j](0,0) = d16 * d19 * d22 + d18 * (d20 * d23 + d21 * d24);
|
||||
Ke[i][j](0,1) = d17 * d19 * d23 + d18 * (d20 * d22);
|
||||
Ke[i][j](0,2) = d17 * d19 * d24 + d18 * (d21 * d22);
|
||||
Ke[i][j](1,0) = d17 * d20 * d22 + d18 * (d19 * d23);
|
||||
Ke[i][j](1,1) = d16 * d20 * d23 + d18 * (d19 * d22 + d21 * d24);
|
||||
Ke[i][j](1,2) = d17 * d20 * d24 + d18 * (d21 * d23);
|
||||
Ke[i][j](2,0) = d17 * d21 * d22 + d18 * (d19 * d24);
|
||||
Ke[i][j](2,1) = d17 * d21 * d23 + d18 * (d20 * d24);
|
||||
Ke[i][j](2,2) = d16 * d21 * d24 + d18 * (d20 * d23 + d19 * d22);
|
||||
Ke[i][j] *= d14;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
/**
|
||||
* Compute Stiffness matrix of tetrahedral element.
|
||||
*
|
||||
*/
|
||||
template<typename real_type, typename vector3_type, typename matrix3x3_type>
|
||||
inline void compute_Ke(
|
||||
vector3_type * B,
|
||||
vector3_type const& D,
|
||||
real_type const & volume,
|
||||
matrix3x3_type Ke[4][4]
|
||||
)
|
||||
{
|
||||
for (unsigned int i=0; i<4; ++i)
|
||||
for (unsigned int j=i; j<4; ++j)
|
||||
compute_Ke_ij( B[i], D, B[j], volume, Ke[i][j] );
|
||||
|
||||
for (unsigned int i=1; i<4; ++i)
|
||||
for (unsigned int j=0; j<i; ++j)
|
||||
Ke[i][j] = trans(Ke[j][i]);
|
||||
}
|
||||
|
||||
/**
|
||||
* Compute i,j sub-block of the element stiffness matrix Ke
|
||||
*
|
||||
* @param Bi Sub-block entries of B matrix of i'th node given in vector-form. Ie. Bi = [bi ci di]^T.
|
||||
* @param D The elasticity matrix in vector form D = [D0 D1 D2]^T.
|
||||
* @param Bj Sub-block entries of B matrix of i'th node given in vector-form. Ie. Bj = [bj cj dj]^T.
|
||||
* @param V The volume of the tetrahedron.
|
||||
* @param Ke_ij Upon return, contains the computed value.
|
||||
*/
|
||||
template<typename real_type, typename vector3_type, typename matrix3x3_type>
|
||||
inline void compute_Ke_ij(
|
||||
vector3_type const& Bi,
|
||||
vector3_type const& D,
|
||||
vector3_type const& Bj,
|
||||
real_type const & volume,
|
||||
matrix3x3_type & Ke_ij
|
||||
)
|
||||
{
|
||||
assert(Ke_ij.size1() == 3 || !"compute_Ke_ij(): Incompatible dimension");
|
||||
assert(Ke_ij.size2() == 3 || !"compute_Ke_ij(): Incompatible dimension");
|
||||
|
||||
//
|
||||
// Calculate Ke_ij = Bi^T D Bj V
|
||||
//
|
||||
// The matrix Bi is given by:
|
||||
//
|
||||
// | bi 0 0 |
|
||||
// | 0 ci 0 |
|
||||
// Bi = | 0 0 di |
|
||||
// | ci bi 0 |
|
||||
// | di 0 bi |
|
||||
// | 0 di ci |
|
||||
//
|
||||
// Similar for Bj. The elasticity matrix D is given by
|
||||
//
|
||||
//
|
||||
// - - - -
|
||||
// | 1-nu nu nu 0 0 0 | | D0 D1 D1 0 0 0 |
|
||||
// young | nu 1-nu nu 0 0 0 | | D1 D0 D1 0 0 0 |
|
||||
// D= ------------- | nu nu 1-nu 0 0 0 | | D1 D1 D0 0 0 0 |
|
||||
// (1+nu)(1-2nu) | 0 0 0 (1-2nu)/2 0 0 | = | 0 0 0 D2 0 0 |
|
||||
// | 0 0 0 0 (1-2nu)/2 0 | | 0 0 0 0 D2 0 |
|
||||
// | 0 0 0 0 0 (1-2nu)/2 | | 0 0 0 0 0 D2 |
|
||||
// - - - -
|
||||
//
|
||||
//
|
||||
// We have something like
|
||||
//
|
||||
// Bi^T D Bj
|
||||
//
|
||||
// These matrices have a particular pattern, especially the B-matrices:
|
||||
//
|
||||
//
|
||||
// |x 0 0|T
|
||||
// |0 x 0|
|
||||
// |0 0 x|
|
||||
// | x 0 0 x 0 x| |x x 0|
|
||||
// | 0 x 0 x x 0| = |0 x x| Bk(i,j)^T = Bk(j,i)
|
||||
// | 0 0 x 0 x x| |x 0 x|
|
||||
//
|
||||
// And the product has the pattern
|
||||
//
|
||||
//
|
||||
// | x x x 0 0 0 | |x 0 0|
|
||||
// | x x x 0 0 0 | |0 x 0|
|
||||
// | x x x 0 0 0 | |0 0 x|
|
||||
// | x 0 0 x 0 x| | 0 0 0 x 0 0 | |x x 0|
|
||||
// | 0 x 0 x x 0| | 0 0 0 0 x 0 | |0 x x|
|
||||
// | 0 0 x 0 x x| | 0 0 0 0 0 x | |x 0 x|
|
||||
//
|
||||
//
|
||||
// If we compute T = Bi^T E, then we get a matrix with pattern
|
||||
//
|
||||
// | x x x x 0 x|
|
||||
// | x x x x x 0|
|
||||
// | x x x 0 x x|
|
||||
//
|
||||
// Now computing Kij = T Bj,
|
||||
//
|
||||
// |x 0 0 |
|
||||
// |0 x 0 |
|
||||
// |0 0 x |
|
||||
// |x x x | | x x x x 0 x| |x x 0 |
|
||||
// |x x x | = | x x x x x 0| |0 x x |
|
||||
// |x x x | | x x x 0 x x| |x 0 x |
|
||||
//
|
||||
//
|
||||
// and exploiting that for isotropic materials we have,
|
||||
//
|
||||
// D(0,0) = D(1,1) = D(2,2) = D0
|
||||
// D(0,1) = D(1,2) = D(0,2) = D(1,0) = D(2,1) = D(2,0) = D1
|
||||
// D(3,3) = D(4,4) = D(5,5) = D2
|
||||
//
|
||||
// Finally we do not need to store the B-matrices, observe that since
|
||||
//
|
||||
// b(0) = invC(1,0); b(1) = invC(1,1); b(2) = invC(1,2); b(3) = invC(1,3);
|
||||
// c(0) = invC(2,0); c(1) = invC(2,1); c(2) = invC(2,2); c(3) = invC(2,3);
|
||||
// d(0) = invC(3,0); d(1) = invC(3,1); d(2) = invC(3,2); d(3) = invC(3,3);
|
||||
//
|
||||
// and
|
||||
//
|
||||
// | bi 0 0 |
|
||||
// | 0 ci 0 |
|
||||
// Bi = | 0 0 di |
|
||||
// | ci bi 0 |
|
||||
// | 0 di ci |
|
||||
// | di 0 bi |
|
||||
//
|
||||
// therefore
|
||||
//
|
||||
// | invC(1,i) 0 0 |
|
||||
// | 0 invC(2,i) 0 |
|
||||
// Bi = | 0 0 invC(3,i) |
|
||||
// | invC(2,i) invC(1,i) 0 |
|
||||
// | 0 invC(3,i) invC(2,i) |
|
||||
// | invC(3,i) 0 invC(1,i) |
|
||||
//
|
||||
// Putting it all together and using matlab, we have,
|
||||
//
|
||||
//syms bi ci di real
|
||||
//Bi = [ bi 0 0 ; 0 ci 0 ; 0 0 di ; ci bi 0 ; 0 di ci ; di 0 bi ]
|
||||
//syms bj cj dj real
|
||||
//Bj = [ bj 0 0 ; 0 cj 0 ; 0 0 dj ; cj bj 0 ; 0 dj cj ; dj 0 bj ]
|
||||
//syms G H J real
|
||||
//D = [ G H H 0 0 0; H G H 0 0 0; H H G 0 0 0; 0 0 0 J 0 0; 0 0 0 0 J 0; 0 0 0 0 0 J]
|
||||
//T = Bi'*E
|
||||
//Kij = T*Bj
|
||||
//syms invC1i invC2i invC3i real;
|
||||
//Kij = subs(Kij,bi,invC1i);
|
||||
//Kij = subs(Kij,ci,invC2i);
|
||||
//Kij = subs(Kij,di,invC3i);
|
||||
//syms invC1j invC2j invC3j real;
|
||||
//Kij = subs(Kij,bj,invC1j);
|
||||
//Kij = subs(Kij,cj,invC2j);
|
||||
//Kij = subs(Kij,dj,invC3j);
|
||||
//Kij = collect(Kij,G)
|
||||
//Kij = collect(Kij,H)
|
||||
//Kij = collect(Kij,J)
|
||||
//
|
||||
// which yields
|
||||
//
|
||||
// - -
|
||||
// | D0 bi bj + D2 (ci cj + di dj) D1 bi cj + D2 ci bj D1 bi dj + D2 di bj |
|
||||
// Ke_ij = V | D1 ci bj + D2 bi cj D0 ci cj + D2 (bi bj + di dj) D1 ci dj + D2 di cj |
|
||||
// | D1 di bj + D2 bi dj D1 di cj + D2 ci dj D0 di dj + D2 (bi bj + ci cj) |
|
||||
// - -
|
||||
real_type bi = Bi(0);
|
||||
real_type ci = Bi(1);
|
||||
real_type di = Bi(2);
|
||||
|
||||
real_type bj = Bj(0);
|
||||
real_type cj = Bj(1);
|
||||
real_type dj = Bj(2);
|
||||
|
||||
real_type D0 = D(0)*volume;
|
||||
real_type D1 = D(1)*volume;
|
||||
real_type D2 = D(2)*volume;
|
||||
|
||||
Ke_ij(0,0) = D0 * bi * bj + D2 * (ci * cj + di * dj);
|
||||
Ke_ij(0,1) = D1 * bi * cj + D2 * (ci * bj);
|
||||
Ke_ij(0,2) = D1 * bi * dj + D2 * (di * bj);
|
||||
|
||||
Ke_ij(1,0) = D1 * ci * bj + D2 * (bi * cj);
|
||||
Ke_ij(1,1) = D0 * ci * cj + D2 * (bi * bj + di * dj);
|
||||
Ke_ij(1,2) = D1 * ci * dj + D2 * (di * cj);
|
||||
|
||||
Ke_ij(2,0) = D1 * di * bj + D2 * (bi * dj);
|
||||
Ke_ij(2,1) = D1 * di * cj + D2 * (ci * dj);
|
||||
Ke_ij(2,2) = D0 * di * dj + D2 * (ci * cj + bi * bj);
|
||||
}
|
||||
|
||||
} // namespace detail
|
||||
} // namespace fem
|
||||
} // namespace OpenTissue
|
||||
|
||||
//OPENTISSUE_DYNAMICS_FEM_FEM_COMPUTE_KE_H
|
||||
#endif
|
||||
@@ -0,0 +1,175 @@
|
||||
#ifndef OPENTISSUE_DYNAMICS_FEM_FEM_COMPUTE_MASS_H
|
||||
#define OPENTISSUE_DYNAMICS_FEM_FEM_COMPUTE_MASS_H
|
||||
//
|
||||
// OpenTissue Template Library
|
||||
// - A generic toolbox for physics-based modeling and simulation.
|
||||
// Copyright (C) 2008 Department of Computer Science, University of Copenhagen.
|
||||
//
|
||||
// OTTL is licensed under zlib: http://opensource.org/licenses/zlib-license.php
|
||||
//
|
||||
#include <OpenTissue/configuration.h>
|
||||
|
||||
#include <limits>
|
||||
|
||||
namespace OpenTissue
|
||||
{
|
||||
namespace fem
|
||||
{
|
||||
namespace detail
|
||||
{
|
||||
|
||||
/**
|
||||
* 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.
|
||||
*
|
||||
* @param mesh
|
||||
*/
|
||||
template<typename fem_mesh>
|
||||
inline void compute_mass(fem_mesh & mesh)
|
||||
{
|
||||
typedef typename fem_mesh::real_type real_type;
|
||||
|
||||
//---
|
||||
//--- A mass matrix is a discrete representation of a continuous mass distribution.
|
||||
//--- To compute our mass matrix for a tetrahedral element with linear shape
|
||||
//--- functions we need the formula (pp. 266 in Cook)
|
||||
//---
|
||||
//--- a!b!c!d!
|
||||
//--- \int_V N_1^a N_2^b N_3^c N_4^d dV = 6V -------------------------- (**)
|
||||
//--- (3 + a + b +c + d)!
|
||||
//---
|
||||
//--- A consistent element mass matrix (pp. 376 Cook) is defined as
|
||||
//---
|
||||
//--- m = \int_V \rho N^T N dV (***)
|
||||
//---
|
||||
//--- This equation can be derived from work balance, the details of which is unimportant
|
||||
//--- here (look Cook pp. 375-376 for details).
|
||||
//--- Assumping \rho is constant over each tetrahedral element and using the linear shape
|
||||
//--- functions the above definition (***) results in
|
||||
//---
|
||||
//--- |N_1|
|
||||
//--- m = \rho \int_V |N_2| |N_1 N_2 N_3 N_4| dV
|
||||
//--- |N_3|
|
||||
//--- |N_4|
|
||||
//---
|
||||
//--- |(N_1 N_1) (N_1 N_2) (N_1 N_3) (N_1 N_4)|
|
||||
//--- m = \rho \int_V |(N_2 N_1) (N_2 N_2) (N_2 N_3) (N_2 N_4)| dV
|
||||
//--- |(N_3 N_1) (N_3 N_2) (N_3 N_3) (N_3 N_4)|
|
||||
//--- |(N_4 N_1) (N_4 N_2) (N_4 N_3) (N_4 N_4)|
|
||||
//---
|
||||
//--- by (**)
|
||||
//---
|
||||
//--- | 2 1 1 1|
|
||||
//--- m = \rho V/20 | 1 2 1 1| (****)
|
||||
//--- | 1 1 2 1|
|
||||
//--- | 1 1 1 2|
|
||||
//---
|
||||
//--- V
|
||||
//--- m_ij = \rho --- (1+delta_ij)
|
||||
//--- 20
|
||||
//---
|
||||
//--- in 3D this means that for the tetrahedral element
|
||||
//---
|
||||
//--- | 2 2 2 1 1 1 1 1 1 1 1 1 |
|
||||
//--- | 2 2 2 1 1 1 1 1 1 1 1 1 |
|
||||
//--- | 2 2 2 1 1 1 1 1 1 1 1 1 |
|
||||
//--- | |
|
||||
//--- | 1 1 1 2 2 2 1 1 1 1 1 1 |
|
||||
//--- | 1 1 1 2 2 2 1 1 1 1 1 1 |
|
||||
//--- V | 1 1 1 2 2 2 1 1 1 1 1 1 |
|
||||
//--- Me = \rho --- | |
|
||||
//--- 20 | 1 1 1 1 1 1 2 2 2 1 1 1 |
|
||||
//--- | 1 1 1 1 1 1 2 2 2 1 1 1 |
|
||||
//--- | 1 1 1 1 1 1 2 2 2 1 1 1 |
|
||||
//--- | |
|
||||
//--- | 1 1 1 1 1 1 1 1 1 2 2 2 |
|
||||
//--- | 1 1 1 1 1 1 1 1 1 2 2 2 |
|
||||
//--- | 1 1 1 1 1 1 1 1 1 2 2 2 |
|
||||
//---
|
||||
//--- Notice that in order to obtain the global/system mass matrix an assembly similar to the
|
||||
//--- stiffnees matrix assembly must be carried out. Further, the global M matrix will
|
||||
//--- have the same sub-block pattern as the global K matrix.
|
||||
//---
|
||||
//--- A consistent mass matrix is often not used in computer graphics. Instead and
|
||||
//--- ad-hoc approach named ``lumped'' mass matrix is applied.
|
||||
//--- The lumped mass matrix is obtained by placing particle masses at the nodes.
|
||||
//--- This corresponds to shifting all the masses in the rows of (****) onto the
|
||||
//--- diagonal. In 3D this yields the element mass matrix
|
||||
//---
|
||||
//--- | 1 0 0 0 0 0 0 0 0 0 0 0 |
|
||||
//--- | 0 1 0 0 0 0 0 0 0 0 0 0 |
|
||||
//--- | 0 0 1 0 0 0 0 0 0 0 0 0 |
|
||||
//--- | |
|
||||
//--- | 0 0 0 1 0 0 0 0 0 0 0 0 |
|
||||
//--- | 0 0 0 0 1 0 0 0 0 0 0 0 |
|
||||
//--- V | 0 0 0 0 0 1 0 0 0 0 0 0 |
|
||||
//--- Me = \rho --- | |
|
||||
//--- 4 | 0 0 0 0 0 0 1 0 0 0 0 0 |
|
||||
//--- | 0 0 0 0 0 0 0 1 0 0 0 0 |
|
||||
//--- | 0 0 0 0 0 0 0 0 1 0 0 0 |
|
||||
//--- | |
|
||||
//--- | 0 0 0 0 0 0 0 0 0 1 0 0 |
|
||||
//--- | 0 0 0 0 0 0 0 0 0 0 1 0 |
|
||||
//--- | 0 0 0 0 0 0 0 0 0 0 0 1 |
|
||||
//---
|
||||
//--- Thus a lumped mass matrix is diagonal whereas a consistent mass matrix
|
||||
//--- is not. Observe that the global mass matrix would also diagonal and the
|
||||
//--- assembly is simplified to an iteration over all tetrahedra, while
|
||||
//--- incementing the nodal mass by one fourth of the tetrahedral mass.
|
||||
//---
|
||||
//--- for each node n
|
||||
//--- mass(n) = 0
|
||||
//--- next n
|
||||
//--- for each tetrahedron e
|
||||
//--- mass(n_i) += \rho_e Ve / 4
|
||||
//--- mass(n_j) += \rho_e Ve / 4
|
||||
//--- mass(n_k) += \rho_e Ve / 4
|
||||
//--- mass(n_m) += \rho_e Ve / 4
|
||||
//--- next e
|
||||
//---
|
||||
//--- where n_i,n_j,n_k and n_m are the four nodes of the e'th tetrahedron.
|
||||
//---
|
||||
//--- The advantage of lumping is less storage and higher performace. On the downside
|
||||
//--- lumping introduces a discontinouty in the displacement field.
|
||||
//---
|
||||
//--- Obrien.shen state that the errors in lumping is negligeble for small-size course
|
||||
//--- meshes used in computer graphics. However, for finer meshes the errors becomes
|
||||
//--- noticeable.
|
||||
//---
|
||||
//--- There do exist other approaches for computing mass matrices, even mehtods which
|
||||
//--- combine other methods. We refer the interested reader to Cook for more details. Here
|
||||
//--- we have limited our selfes to the two most common methods.
|
||||
//---
|
||||
//--- It is worthwhile to notice that under the reasonable assumptions that V and \rho are
|
||||
//--- positive for all elements both the element mass matrices and the global mass matrices
|
||||
//--- are symmetric positive definite matrices.
|
||||
|
||||
for (int n=0;n<mesh.m_nodes.size();n++)
|
||||
{
|
||||
|
||||
if(mesh.m_nodes[n].m_fixed)
|
||||
mesh.m_nodes[n].m_mass = std::numeric_limits<real_type>::max();
|
||||
else
|
||||
mesh.m_nodes[n].m_mass = 0;
|
||||
}
|
||||
|
||||
for (int t=0;t<mesh.m_tetrahedra.size();t++)
|
||||
{
|
||||
|
||||
real_type amount = mesh.m_tetrahedra[t].m_density * mesh.m_tetrahedra[t].m_V *.25;
|
||||
mesh.m_tetrahedra[t].i()->m_mass += amount;
|
||||
mesh.m_tetrahedra[t].j()->m_mass += amount;
|
||||
mesh.m_tetrahedra[t].k()->m_mass += amount;
|
||||
mesh.m_tetrahedra[t].m()->m_mass += amount;
|
||||
}
|
||||
}
|
||||
|
||||
} // namespace detail
|
||||
} // namespace fem
|
||||
} // namespace OpenTissue
|
||||
|
||||
//OPENTISSUE_DYNAMICS_FEM_FEM_COMPUTE_MASS_H
|
||||
#endif
|
||||
@@ -0,0 +1,93 @@
|
||||
#ifndef OPENTISSUE_DYNAMICS_FEM_FEM_COMPUTE_VOLUME_H
|
||||
#define OPENTISSUE_DYNAMICS_FEM_FEM_COMPUTE_VOLUME_H
|
||||
//
|
||||
// OpenTissue Template Library
|
||||
// - A generic toolbox for physics-based modeling and simulation.
|
||||
// Copyright (C) 2008 Department of Computer Science, University of Copenhagen.
|
||||
//
|
||||
// OTTL is licensed under zlib: http://opensource.org/licenses/zlib-license.php
|
||||
//
|
||||
#include <OpenTissue/configuration.h>
|
||||
|
||||
namespace OpenTissue
|
||||
{
|
||||
namespace fem
|
||||
{
|
||||
namespace detail
|
||||
{
|
||||
/**
|
||||
* Compute volume of tetrahedron.
|
||||
*
|
||||
* @param e10 edge vector from p0 to p1, i.e. e10 = p1-p0.
|
||||
* @param e20 edge vector from p0 to p2, i.e. e20 = p2-p0.
|
||||
* @param e30 edge vector from p0 to p3, i.e. e30 = p3-p0.
|
||||
*
|
||||
* @return The signed volume of the tetrahedra with nodes p0,p1,p2 and p3.
|
||||
*/
|
||||
template<typename vector3_type>
|
||||
inline typename vector3_type::value_type compute_volume(vector3_type const& e10, vector3_type const& e20, vector3_type const& e30)
|
||||
{
|
||||
typedef typename vector3_type::value_type real_type;
|
||||
// Calculate the volume given by:
|
||||
// | 1 1 1 1 |
|
||||
// | e1x e2x e3x | | p0x p1x p2x p3x |
|
||||
// 6V = | e1 \cdot (e2 x e3) | = | e1y e2y e3y | = 6V | p0y p1y p2y p3y |
|
||||
// | e1z e2z e3z | | p0z p1z p2z p3z |
|
||||
//
|
||||
real_type sixV = e10*(e20 % e30);
|
||||
return sixV / 6.0;
|
||||
}
|
||||
|
||||
/**
|
||||
* Compute Tetrahedron Volume
|
||||
*
|
||||
* @param p0
|
||||
* @param p1
|
||||
* @param p2
|
||||
* @param p3
|
||||
*
|
||||
* @return The signed volume of the tetrahedra with nodes p0,p1,p2 and p3.
|
||||
*/
|
||||
template<typename vector3_type>
|
||||
inline typename vector3_type::value_type compute_volume(vector3_type const & p0,vector3_type const & p1,vector3_type const & p2,vector3_type const & p3)
|
||||
{
|
||||
typedef typename vector3_type::value_type real_type;
|
||||
|
||||
real_type d = p0(0);
|
||||
real_type d4 = p0(1);
|
||||
real_type d8 = p0(2);
|
||||
real_type d1 = p1(0) - d;
|
||||
real_type d5 = p1(1) - d4;
|
||||
real_type d9 = p1(2) - d8;
|
||||
real_type d2 = p2(0) - d;
|
||||
real_type d6 = p2(1) - d4;
|
||||
real_type d10 = p2(2) - d8;
|
||||
real_type d3 = p3(0) - d;
|
||||
real_type d7 = p3(1) - d4;
|
||||
real_type d11 = p3(2) - d8;
|
||||
real_type d12 = (d1 * d6 * d11 + d2 * d7 * d9 + d3 * d5 * d10) - d1 * d7 * d10 - d2 * d5 * d11 - d3 * d6 * d9;
|
||||
return d12 / 6.0;
|
||||
/*
|
||||
** From Zienkiewicz & Taylor p.637
|
||||
** V = 1/6*det( 1 x0 y0 z0; 1 x1 y1 z1; 1 x2 y2 z2; 1 x3 y3 z3)
|
||||
** where x0 = n0->i; y0 = n0(1); ...
|
||||
** Calculated by Mathematica ;-)
|
||||
*/
|
||||
/*
|
||||
return (p0(2)*p1(1)*p2(0) - p0(1)*p1(2)*p2(0) - p0(2)*p1(0)*p2(1)
|
||||
+ p0(0)*p1(2)*p2(1) + p0(1)*p1(0)*p2(2) - p0(0)*p1(1)*p2(2)
|
||||
- p0(2)*p1(1)*p3(0) + p0(1)*p1(2)*p3(0) + p0(2)*p2(1)*p3(0)
|
||||
- p1(2)*p2(1)*p3(0) - p0(1)*p2(2)*p3(0) + p1(1)*p2(2)*p3(0)
|
||||
+ p0(2)*p1(0)*p3(1) - p0(0)*p1(2)*p3(1) - p0(2)*p2(0)*p3(1)
|
||||
+ p1(2)*p2(0)*p3(1) + p0(0)*p2(2)*p3(1) - p1(0)*p2(2)*p3(1)
|
||||
- p0(1)*p1(0)*p3(2) + p0(0)*p1(1)*p3(2) + p0(1)*p2(0)*p3(2)
|
||||
- p1(1)*p2(0)*p3(2) - p0(0)*p2(1)*p3(2) + p1(0)*p2(1)*p3(2))/6.;
|
||||
*/
|
||||
}
|
||||
|
||||
} // namespace detail
|
||||
} // namespace fem
|
||||
} // namespace OpenTissue
|
||||
|
||||
//OPENTISSUE_DYNAMICS_FEM_FEM_COMPUTE_VOLUME_H
|
||||
#endif
|
||||
@@ -0,0 +1,141 @@
|
||||
#ifndef OPENTISSUE_DYNAMICS_FEM_FEM_CONJUGATE_GRADIENTS_H
|
||||
#define OPENTISSUE_DYNAMICS_FEM_FEM_CONJUGATE_GRADIENTS_H
|
||||
//
|
||||
// OpenTissue Template Library
|
||||
// - A generic toolbox for physics-based modeling and simulation.
|
||||
// Copyright (C) 2008 Department of Computer Science, University of Copenhagen.
|
||||
//
|
||||
// OTTL is licensed under zlib: http://opensource.org/licenses/zlib-license.php
|
||||
//
|
||||
#include <OpenTissue/configuration.h>
|
||||
|
||||
#define n_i (&mesh.m_nodes[n])
|
||||
#define n_j (&mesh.m_nodes[j])
|
||||
|
||||
|
||||
namespace OpenTissue
|
||||
{
|
||||
namespace fem
|
||||
{
|
||||
namespace detail
|
||||
{
|
||||
/**
|
||||
* Conjugate Gradient Solver.
|
||||
* This solver has been hard-wired for a WarpT4Mesh to solve the equation
|
||||
*
|
||||
* A v = b
|
||||
*
|
||||
*/
|
||||
template < typename fem_mesh >
|
||||
inline void conjugate_gradients(
|
||||
fem_mesh & mesh
|
||||
, unsigned int min_iterations
|
||||
, unsigned int max_iterations
|
||||
)
|
||||
{
|
||||
using std::fabs;
|
||||
|
||||
typedef typename fem_mesh::real_type real_type;
|
||||
typedef typename fem_mesh::vector3_type vector3_type;
|
||||
typedef typename fem_mesh::matrix3x3_type matrix3x3_type;
|
||||
typedef typename fem_mesh::node_type::matrix_iterator matrix_iterator;
|
||||
|
||||
|
||||
real_type tiny = 1e-010; // TODO: Should be user controllable
|
||||
real_type tolerence = 0.001; // TODO: Should be user controllable
|
||||
|
||||
//--- r = b - A v
|
||||
//--- p = r
|
||||
for(int n=0;n<mesh.m_nodes.size();n++)
|
||||
{
|
||||
|
||||
if(n_i->m_fixed)
|
||||
continue;
|
||||
|
||||
n_i->m_residual = n_i->m_b;
|
||||
|
||||
matrix_iterator Abegin = n_i->Abegin();
|
||||
matrix_iterator Aend = n_i->Aend();
|
||||
for (matrix_iterator A = Abegin; A != Aend;++A)
|
||||
{
|
||||
unsigned int j = A->first;
|
||||
matrix3x3_type & A_ij = A->second;
|
||||
vector3_type & v_j = n_j->m_velocity;
|
||||
|
||||
n_i->m_residual -= A_ij * v_j;
|
||||
}
|
||||
n_i->m_prev = n_i->m_residual;
|
||||
}
|
||||
|
||||
for(unsigned int iteration = 0; iteration < max_iterations; ++iteration)
|
||||
{
|
||||
real_type d = 0.0;
|
||||
real_type d2 = 0.0;
|
||||
|
||||
//--- u = A p
|
||||
//--- d = r*r
|
||||
//--- d2 = p*u
|
||||
|
||||
for(int n=0;n<mesh.m_nodes.size();n++)
|
||||
{
|
||||
if(n_i->m_fixed)
|
||||
continue;
|
||||
|
||||
n_i->m_update.clear();
|
||||
|
||||
matrix_iterator Abegin = n_i->Abegin();
|
||||
matrix_iterator Aend = n_i->Aend();
|
||||
for (matrix_iterator A = Abegin; A != Aend;++A)
|
||||
{
|
||||
|
||||
unsigned int j = A->first;
|
||||
matrix3x3_type & A_ij = A->second;
|
||||
|
||||
n_i->m_update += A_ij * n_j->m_prev;
|
||||
|
||||
}
|
||||
d += n_i->m_residual * n_i->m_residual;
|
||||
d2 += n_i->m_prev * n_i->m_update;
|
||||
}
|
||||
|
||||
if(fabs(d2) < tiny)
|
||||
d2 = tiny;
|
||||
|
||||
real_type d3 = d / d2;
|
||||
real_type d1 = 0.0;
|
||||
//--- v += p * d3
|
||||
//--- r -= u * d3
|
||||
//--- d1 = r*r
|
||||
|
||||
for(int n=0;n<mesh.m_nodes.size();n++)
|
||||
{
|
||||
if(n_i->m_fixed)
|
||||
continue;
|
||||
n_i->m_velocity += n_i->m_prev * d3;
|
||||
n_i->m_residual -= n_i->m_update * d3;
|
||||
d1 += n_i->m_residual * n_i->m_residual;
|
||||
}
|
||||
if(iteration >= min_iterations && d1 < tolerence)
|
||||
break;
|
||||
|
||||
if(fabs(d) < tiny)
|
||||
d = tiny;
|
||||
|
||||
real_type d4 = d1 / d;
|
||||
//--- p = r + d4 * p
|
||||
for(int n=0;n<mesh.m_nodes.size();n++)
|
||||
{
|
||||
|
||||
if(n_i->m_fixed)
|
||||
continue;
|
||||
n_i->m_prev = n_i->m_residual + n_i->m_prev * d4;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
} // namespace detail
|
||||
} // namespace fem
|
||||
} // namespace OpenTissue
|
||||
|
||||
//OPENTISSUE_DYNAMICS_FEM_FEM_CONJUGATE_GRADIENTS_H
|
||||
#endif
|
||||
@@ -0,0 +1,104 @@
|
||||
#ifndef OPENTISSUE_DYNAMICS_FEM_FEM_DYNAMICS_ASSEMBLY_H
|
||||
#define OPENTISSUE_DYNAMICS_FEM_FEM_DYNAMICS_ASSEMBLY_H
|
||||
//
|
||||
// OpenTissue Template Library
|
||||
// - A generic toolbox for physics-based modeling and simulation.
|
||||
// Copyright (C) 2008 Department of Computer Science, University of Copenhagen.
|
||||
//
|
||||
// OTTL is licensed under zlib: http://opensource.org/licenses/zlib-license.php
|
||||
//
|
||||
#include <OpenTissue/configuration.h>
|
||||
|
||||
#define n_i (&mesh.m_nodes[n])
|
||||
#define n_j (&mesh.m_nodes[j])
|
||||
|
||||
|
||||
namespace OpenTissue
|
||||
{
|
||||
namespace fem
|
||||
{
|
||||
namespace detail
|
||||
{
|
||||
/**
|
||||
* Setup dynamic equation.
|
||||
* This is the equation of motion given as:
|
||||
*
|
||||
*
|
||||
* A v^{i+1} = b
|
||||
*
|
||||
*
|
||||
* where
|
||||
*
|
||||
|
||||
* A = M + \delta t C + \delta t^2 K
|
||||
* b = M v^i - \delta 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.
|
||||
*
|
||||
* @param dt The time step, \delta t, which is about to be taken.
|
||||
* @param mass_damping Coefficient for mass damping in the Raleigh damping equation.
|
||||
* The coefficient \alpha in C = \alpha M + \beta K. In this implementation \beta = 0.
|
||||
*
|
||||
*
|
||||
*
|
||||
*
|
||||
*/
|
||||
template < typename fem_mesh, typename real_type >
|
||||
inline void dynamics_assembly(
|
||||
fem_mesh & mesh,
|
||||
real_type const & mass_damping,
|
||||
real_type const & dt
|
||||
)
|
||||
{
|
||||
typedef typename fem_mesh::vector3_type vector3_type;
|
||||
typedef typename fem_mesh::matrix3x3_type matrix3x3_type;
|
||||
typedef typename fem_mesh::node_type::matrix_iterator matrix_iterator;
|
||||
|
||||
for (int n = 0;n<mesh.m_nodes.size();n++)
|
||||
{
|
||||
fem_mesh::node_type& ni = mesh.m_nodes[n];
|
||||
if (ni.m_mass)
|
||||
{
|
||||
unsigned int i = n_i->idx();
|
||||
vector3_type & b_i = n_i->m_b;
|
||||
real_type & m_i = n_i->m_mass;
|
||||
|
||||
b_i.clear();
|
||||
|
||||
matrix_iterator Kbegin = n_i->Kbegin();
|
||||
matrix_iterator Kend = n_i->Kend();
|
||||
for (matrix_iterator K = Kbegin; K != Kend;++K)
|
||||
{
|
||||
unsigned int j = K->first;
|
||||
matrix3x3_type & K_ij = K->second;
|
||||
vector3_type & x_j = n_j->m_coord;
|
||||
matrix3x3_type & A_ij = n_i->A(j);
|
||||
|
||||
A_ij = K_ij * (dt*dt);
|
||||
b_i -= K_ij * x_j;
|
||||
if (i == j)
|
||||
{
|
||||
real_type c_i = mass_damping*m_i;
|
||||
real_type tmp = m_i + dt*c_i;
|
||||
A_ij(0,0) += tmp; A_ij(1,1) += tmp; A_ij(2,2) += tmp;
|
||||
}
|
||||
}
|
||||
b_i -= n_i->m_f0;
|
||||
b_i += n_i->m_f_external;
|
||||
b_i *= dt;
|
||||
b_i += n_i->m_velocity * m_i;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
} // namespace detail
|
||||
} // namespace fem
|
||||
} // namespace OpenTissue
|
||||
|
||||
//OPENTISSUE_DYNAMICS_FEM_FEM_DYNAMICS_ASSEMBLY_H
|
||||
#endif
|
||||
@@ -0,0 +1,84 @@
|
||||
#ifndef OPENTISSUE_DYNAMICS_FEM_FEM_INIT_H
|
||||
#define OPENTISSUE_DYNAMICS_FEM_FEM_INIT_H
|
||||
//
|
||||
// OpenTissue Template Library
|
||||
// - A generic toolbox for physics-based modeling and simulation.
|
||||
// Copyright (C) 2008 Department of Computer Science, University of Copenhagen.
|
||||
//
|
||||
// OTTL is licensed under zlib: http://opensource.org/licenses/zlib-license.php
|
||||
//
|
||||
#include <OpenTissue/configuration.h>
|
||||
|
||||
#include <OpenTissue/dynamics/fem/fem_uniform_young.h>
|
||||
#include <OpenTissue/dynamics/fem/fem_uniform_poisson.h>
|
||||
#include <OpenTissue/dynamics/fem/fem_uniform_density.h>
|
||||
#include <OpenTissue/dynamics/fem/fem_compute_mass.h>
|
||||
#include <OpenTissue/dynamics/fem/fem_initialize_stiffness_elements.h>
|
||||
#include <OpenTissue/dynamics/fem/fem_clear_stiffness_assembly.h>
|
||||
#include <OpenTissue/dynamics/fem/fem_initialize_plastic.h>
|
||||
|
||||
namespace OpenTissue
|
||||
{
|
||||
namespace fem
|
||||
{
|
||||
/**
|
||||
* Initialize Data Structres.
|
||||
* Should be invoked at start-up before any calls to animate_warping.
|
||||
*
|
||||
* This function assumes that mesh geometry have been set up prior
|
||||
* to invokation (i.e. both world coord and orignial coord are initlized
|
||||
* and the same).
|
||||
*
|
||||
* This method initialized the model with the same material parameters. End
|
||||
* user could write his or her own initilization method replacing
|
||||
* the steps where young modules, poisson ratio and densities are
|
||||
* assigned to the model.
|
||||
*
|
||||
* @param mesh
|
||||
* @param young Young Moduls, a value of 500000 seems work quite well.
|
||||
* @param poisson Poisson ration, a value of 0.33 seems work quite well.
|
||||
* @param density Mass density, a value of 1000 seems to work quite well.
|
||||
* @param c_yield Plastic yield.
|
||||
* @param c_creep Plastic creep.
|
||||
* @param c_max Plastic max.
|
||||
*/
|
||||
template < typename fem_mesh,typename real_type >
|
||||
inline void init(
|
||||
fem_mesh & mesh,
|
||||
real_type const & young,
|
||||
real_type const & poisson,
|
||||
real_type const & density,
|
||||
real_type const & c_yield,
|
||||
real_type const & c_creep,
|
||||
real_type const & c_max
|
||||
)
|
||||
{
|
||||
assert(poisson>=0 || !"uniform_poisson() : Poisson ratio must be non-negative");
|
||||
assert(poisson<=0.5 || !"uniform_poisson() : Poisson ratio must not be larger than a half");
|
||||
assert(density>0 || !"uniform_density(): density must be positive");
|
||||
|
||||
//--- Assign material parameters
|
||||
for (int t=0;t<mesh.m_tetrahedra.size();t++)
|
||||
{
|
||||
mesh.m_tetrahedra[t].m_young = young;
|
||||
mesh.m_tetrahedra[t].m_poisson = poisson;
|
||||
mesh.m_tetrahedra[t].m_density = density;
|
||||
detail::initialize_stiffness_elements_single(&mesh.m_tetrahedra[t]);
|
||||
}
|
||||
//--- Compute stiffness and mass matrices
|
||||
for (int m=0;m<mesh.m_nodes.size();m++)
|
||||
{
|
||||
detail::clear_stiffness_assembly_single(&mesh.m_nodes[m]);
|
||||
}
|
||||
detail::compute_mass(mesh);
|
||||
for (int t=0;t<mesh.m_tetrahedra.size();t++)
|
||||
{
|
||||
detail::initialize_plastic_single(&mesh.m_tetrahedra[t],c_yield,c_creep,c_max);
|
||||
}
|
||||
}
|
||||
|
||||
} // namespace fem
|
||||
} // namespace OpenTissue
|
||||
|
||||
//OPENTISSUE_DYNAMICS_FEM_FEM_INIT_H
|
||||
#endif
|
||||
@@ -0,0 +1,60 @@
|
||||
#ifndef OPENTISSUE_DYNAMICS_FEM_FEM_INITIALIZE_PLASTIC_H
|
||||
#define OPENTISSUE_DYNAMICS_FEM_FEM_INITIALIZE_PLASTIC_H
|
||||
//
|
||||
// OpenTissue Template Library
|
||||
// - A generic toolbox for physics-based modeling and simulation.
|
||||
// Copyright (C) 2008 Department of Computer Science, University of Copenhagen.
|
||||
//
|
||||
// OTTL is licensed under zlib: http://opensource.org/licenses/zlib-license.php
|
||||
//
|
||||
#include <OpenTissue/configuration.h>
|
||||
|
||||
#include <OpenTissue/dynamics/fem/fem_compute_b.h>
|
||||
#include <OpenTissue/dynamics/fem/fem_compute_isotropic_elasticity.h>
|
||||
|
||||
namespace OpenTissue
|
||||
{
|
||||
namespace fem
|
||||
{
|
||||
namespace detail
|
||||
{
|
||||
/**
|
||||
* Initialize plasticity material parameters.
|
||||
*
|
||||
* Note: Initialize_stiffness_elements must have been invoked prior
|
||||
* to calling this function.
|
||||
*
|
||||
* @param begin Iterator to first tetrahedron.
|
||||
* @param end Iterator to one past last tetrahedron.
|
||||
* @param c_yield Plastic yield.
|
||||
* @param c_creep Plastic creep.
|
||||
* @param c_max Plastic max.
|
||||
*/
|
||||
template < typename tetrahedron_type,typename real_type >
|
||||
inline void initialize_plastic_single(
|
||||
tetrahedron_type* T
|
||||
, real_type const & c_yield
|
||||
, real_type const & c_creep
|
||||
, real_type const & c_max
|
||||
)
|
||||
{
|
||||
assert(c_yield>=0 || !"initialize_plastic(): yield must be non-negative");
|
||||
assert(c_creep>=0 || !"initialize_plastic(): creep must be non-negative");
|
||||
assert(c_max>=0 || !"initialize_plastic(): max must be non-negative");
|
||||
|
||||
T->m_yield = c_yield;
|
||||
T->m_creep = c_creep;
|
||||
T->m_max = c_max;
|
||||
for(int i=0;i<6;++i)
|
||||
T->m_plastic[i] = 0;
|
||||
|
||||
compute_B(T->m_e10, T->m_e20, T->m_e30, T->m_V, T->m_B);
|
||||
compute_isotropic_elasticity_vector (T->m_young, T->m_poisson, T->m_D);
|
||||
}
|
||||
|
||||
} // namespace detail
|
||||
} // namespace fem
|
||||
} // namespace OpenTissue
|
||||
|
||||
//OPENTISSUE_DYNAMICS_FEM_FEM_INITIALIZE_PLASTIC_H
|
||||
#endif
|
||||
@@ -0,0 +1,68 @@
|
||||
#ifndef OPENTISSUE_DYNAMICS_FEM_FEM_INITIALIZE_STIFFNESS_ELEMENTS_H
|
||||
#define OPENTISSUE_DYNAMICS_FEM_FEM_INITIALIZE_STIFFNESS_ELEMENTS_H
|
||||
//
|
||||
// OpenTissue Template Library
|
||||
// - A generic toolbox for physics-based modeling and simulation.
|
||||
// Copyright (C) 2008 Department of Computer Science, University of Copenhagen.
|
||||
//
|
||||
// OTTL is licensed under zlib: http://opensource.org/licenses/zlib-license.php
|
||||
//
|
||||
#include <OpenTissue/configuration.h>
|
||||
|
||||
#include <OpenTissue/dynamics/fem/fem_compute_ke.h>
|
||||
#include <OpenTissue/dynamics/fem/fem_compute_volume.h>
|
||||
|
||||
namespace OpenTissue
|
||||
{
|
||||
namespace fem
|
||||
{
|
||||
namespace detail
|
||||
{
|
||||
/**
|
||||
*
|
||||
*
|
||||
* NOTE: Material parameters (young modulus and poisson ratio) must have
|
||||
* been set prior to invoking this function.
|
||||
*
|
||||
* @param begin
|
||||
* @param end
|
||||
*/
|
||||
template < typename tetrahedron_type >
|
||||
inline void initialize_stiffness_elements_single(tetrahedron_type* T)
|
||||
{
|
||||
{
|
||||
T->m_e10 = T->j()->m_model_coord - T->i()->m_model_coord;
|
||||
T->m_e20 = T->k()->m_model_coord - T->i()->m_model_coord;
|
||||
T->m_e30 = T->m()->m_model_coord - T->i()->m_model_coord;
|
||||
|
||||
T->m_V = compute_volume(T->m_e10, T->m_e20, T->m_e30);
|
||||
|
||||
//T->m_V = compute_volume(
|
||||
// T->i()->m_model_coord,
|
||||
// T->j()->m_model_coord,
|
||||
// T->k()->m_model_coord,
|
||||
// T->m()->m_model_coord
|
||||
// );
|
||||
|
||||
assert(T->m_V>0 || !"initialize_stiffness_elements(): Element with negative volume is encountered! Maybe you got a bad mesh?");
|
||||
|
||||
T->m_Re = math::diag(1.0);
|
||||
|
||||
//compute_Ke(T->m_B, T->m_D, T->m_V, T->m_Ke);
|
||||
compute_Ke(
|
||||
T->node(0)->m_model_coord,
|
||||
T->node(1)->m_model_coord,
|
||||
T->node(2)->m_model_coord,
|
||||
T->node(3)->m_model_coord,
|
||||
T->m_young, T->m_poisson,
|
||||
T->m_Ke
|
||||
);
|
||||
}
|
||||
}
|
||||
|
||||
} // namespace detail
|
||||
} // namespace fem
|
||||
} // namespace OpenTissue
|
||||
|
||||
//OPENTISSUE_DYNAMICS_FEM_FEM_INITIALIZE_STIFFNESS_ELEMENTS_H
|
||||
#endif
|
||||
@@ -0,0 +1,40 @@
|
||||
#ifndef OPENTISSUE_DYNAMICS_FEM_FEM_MESH_H
|
||||
#define OPENTISSUE_DYNAMICS_FEM_FEM_MESH_H
|
||||
//
|
||||
// OpenTissue Template Library
|
||||
// - A generic toolbox for physics-based modeling and simulation.
|
||||
// Copyright (C) 2008 Department of Computer Science, University of Copenhagen.
|
||||
//
|
||||
// OTTL is licensed under zlib: http://opensource.org/licenses/zlib-license.php
|
||||
//
|
||||
#include <OpenTissue/configuration.h>
|
||||
|
||||
#include <OpenTissue/core/containers/t4mesh/t4mesh.h>
|
||||
#include <OpenTissue/dynamics/fem/fem_node_traits.h>
|
||||
#include <OpenTissue/dynamics/fem/fem_tetrahedron_traits.h>
|
||||
|
||||
namespace OpenTissue
|
||||
{
|
||||
namespace fem
|
||||
{
|
||||
|
||||
template <typename math_types>
|
||||
class Mesh
|
||||
: public OpenTissue::t4mesh::T4Mesh<
|
||||
math_types
|
||||
, OpenTissue::fem::detail::NodeTraits<math_types>
|
||||
, OpenTissue::fem::detail::TetrahedronTraits<math_types>
|
||||
>
|
||||
{
|
||||
public:
|
||||
|
||||
typedef typename math_types::real_type real_type;
|
||||
typedef typename math_types::vector3_type vector3_type;
|
||||
typedef typename math_types::matrix3x3_type matrix3x3_type;
|
||||
};
|
||||
|
||||
} // namespace fem
|
||||
} // namespace OpenTissue
|
||||
|
||||
//OPENTISSUE_DYNAMICS_FEM_FEM_MESH_H
|
||||
#endif
|
||||
@@ -0,0 +1,86 @@
|
||||
#ifndef OPENTISSUE_DYNAMICS_FEM_FEM_NODE_TRAITS_H
|
||||
#define OPENTISSUE_DYNAMICS_FEM_FEM_NODE_TRAITS_H
|
||||
//
|
||||
// OpenTissue Template Library
|
||||
// - A generic toolbox for physics-based modeling and simulation.
|
||||
// Copyright (C) 2008 Department of Computer Science, University of Copenhagen.
|
||||
//
|
||||
// OTTL is licensed under zlib: http://opensource.org/licenses/zlib-license.php
|
||||
//
|
||||
#include <OpenTissue/configuration.h>
|
||||
|
||||
#include <map>
|
||||
|
||||
namespace OpenTissue
|
||||
{
|
||||
namespace fem
|
||||
{
|
||||
namespace detail
|
||||
{
|
||||
|
||||
template <typename math_types>
|
||||
class NodeTraits
|
||||
{
|
||||
public:
|
||||
|
||||
typedef typename math_types::real_type real_type;
|
||||
typedef typename math_types::vector3_type vector3_type;
|
||||
typedef typename math_types::matrix3x3_type matrix3x3_type;
|
||||
|
||||
typedef typename std::map<int,matrix3x3_type> matrix_container;
|
||||
typedef typename matrix_container::iterator matrix_iterator;
|
||||
|
||||
public:
|
||||
|
||||
NodeTraits(int i,int j)
|
||||
{
|
||||
}
|
||||
matrix_container m_K_row; ///< Currently stored in a map container, key correspond to column and mapped value to 3-by-3 sub block.
|
||||
matrix_container m_A_row;
|
||||
vector3_type m_f0;
|
||||
vector3_type m_b;
|
||||
|
||||
bool m_fixed; ///< If the node is in-moveable this flag is set to true otherwise it is false.
|
||||
|
||||
vector3_type m_model_coord;
|
||||
vector3_type m_coord; ///< World coord
|
||||
|
||||
vector3_type m_velocity;
|
||||
real_type m_mass;
|
||||
vector3_type m_f_external;
|
||||
|
||||
// Needed by the ConjugateGradient method
|
||||
vector3_type m_update;
|
||||
vector3_type m_prev;
|
||||
vector3_type m_residual;
|
||||
|
||||
public:
|
||||
|
||||
matrix_iterator Kbegin()
|
||||
{
|
||||
return m_K_row.begin();
|
||||
}
|
||||
matrix_iterator Kend() { return m_K_row.end(); }
|
||||
matrix_iterator Abegin() { return m_A_row.begin(); }
|
||||
matrix_iterator Aend() { return m_A_row.end(); }
|
||||
|
||||
matrix3x3_type & K(int const & column_idx) {
|
||||
return m_K_row[column_idx];
|
||||
}
|
||||
matrix3x3_type & A(int const & column_idx) { return m_A_row[column_idx]; }
|
||||
|
||||
public:
|
||||
|
||||
NodeTraits()
|
||||
: m_fixed(false)
|
||||
{
|
||||
}
|
||||
|
||||
};
|
||||
|
||||
} // namespace detail
|
||||
} // namespace fem
|
||||
} // namespace OpenTissue
|
||||
|
||||
//OPENTISSUE_DYNAMICS_FEM_FEM_NODE_TRAITS_H
|
||||
#endif
|
||||
@@ -0,0 +1,41 @@
|
||||
#ifndef OPENTISSUE_DYNAMICS_FEM_FEM_POSITION_UPDATE_H
|
||||
#define OPENTISSUE_DYNAMICS_FEM_FEM_POSITION_UPDATE_H
|
||||
//
|
||||
// OpenTissue Template Library
|
||||
// - A generic toolbox for physics-based modeling and simulation.
|
||||
// Copyright (C) 2008 Department of Computer Science, University of Copenhagen.
|
||||
//
|
||||
// OTTL is licensed under zlib: http://opensource.org/licenses/zlib-license.php
|
||||
//
|
||||
#include <OpenTissue/configuration.h>
|
||||
|
||||
namespace OpenTissue
|
||||
{
|
||||
namespace fem
|
||||
{
|
||||
namespace detail
|
||||
{
|
||||
/**
|
||||
* Position Update.
|
||||
*
|
||||
* Note: Velocities must have been computed prior to invoking this function.
|
||||
* @param mesh
|
||||
* @param dt
|
||||
*/
|
||||
template < typename fem_mesh, typename real_type >
|
||||
inline void position_update(fem_mesh & mesh, real_type const & dt)
|
||||
{
|
||||
for (int n=0;n<mesh.m_nodes.size();n++)
|
||||
{
|
||||
if(mesh.m_nodes[n].m_fixed)
|
||||
continue;
|
||||
mesh.m_nodes[n].m_coord += dt * mesh.m_nodes[n].m_velocity;
|
||||
}
|
||||
}
|
||||
|
||||
} // namespace detail
|
||||
} // namespace fem
|
||||
} // namespace OpenTissue
|
||||
|
||||
//OPENTISSUE_DYNAMICS_FEM_FEM_POSITION_UPDATE_H
|
||||
#endif
|
||||
@@ -0,0 +1,33 @@
|
||||
#ifndef OPENTISSUE_DYNAMICS_FEM_FEM_RESET_ORIENTATION_H
|
||||
#define OPENTISSUE_DYNAMICS_FEM_FEM_RESET_ORIENTATION_H
|
||||
//
|
||||
// OpenTissue Template Library
|
||||
// - A generic toolbox for physics-based modeling and simulation.
|
||||
// Copyright (C) 2008 Department of Computer Science, University of Copenhagen.
|
||||
//
|
||||
// OTTL is licensed under zlib: http://opensource.org/licenses/zlib-license.php
|
||||
//
|
||||
#include <OpenTissue/configuration.h>
|
||||
|
||||
namespace OpenTissue
|
||||
{
|
||||
namespace fem
|
||||
{
|
||||
namespace detail
|
||||
{
|
||||
|
||||
/**
|
||||
*
|
||||
*/
|
||||
template<typename tetrahedron_type>
|
||||
inline void reset_orientation_single(tetrahedron_type* T)
|
||||
{
|
||||
T->m_Re = math::diag(1.0);
|
||||
}
|
||||
|
||||
} // namespace detail
|
||||
} // namespace fem
|
||||
} // namespace OpenTissue
|
||||
|
||||
//OPENTISSUE_DYNAMICS_FEM_FEM_RESET_ORIENTATION_H
|
||||
#endif
|
||||
@@ -0,0 +1,131 @@
|
||||
#ifndef OPENTISSUE_DYNAMICS_FEM_FEM_SIMULATE_H
|
||||
#define OPENTISSUE_DYNAMICS_FEM_FEM_SIMULATE_H
|
||||
//
|
||||
// OpenTissue Template Library
|
||||
// - A generic toolbox for physics-based modeling and simulation.
|
||||
// Copyright (C) 2008 Department of Computer Science, University of Copenhagen.
|
||||
//
|
||||
// OTTL is licensed under zlib: http://opensource.org/licenses/zlib-license.php
|
||||
//
|
||||
#include <OpenTissue/configuration.h>
|
||||
|
||||
#include <OpenTissue/dynamics/fem/fem_clear_stiffness_assembly.h>
|
||||
#include <OpenTissue/dynamics/fem/fem_update_orientation.h>
|
||||
#include <OpenTissue/dynamics/fem/fem_reset_orientation.h>
|
||||
#include <OpenTissue/dynamics/fem/fem_stiffness_assembly.h>
|
||||
#include <OpenTissue/dynamics/fem/fem_add_plasticity_force.h>
|
||||
#include <OpenTissue/dynamics/fem/fem_dynamics_assembly.h>
|
||||
#include <OpenTissue/dynamics/fem/fem_conjugate_gradients.h>
|
||||
#include <OpenTissue/dynamics/fem/fem_position_update.h>
|
||||
|
||||
namespace OpenTissue
|
||||
{
|
||||
namespace fem
|
||||
{
|
||||
/**
|
||||
* Simulate.
|
||||
* Note external forces must have been computed prior to invocation.
|
||||
*
|
||||
* The dynamic equation has the following form (where u=x-x_0, and x' is derivative wrt. time)
|
||||
*
|
||||
* M x'' + Cx' + K (x-x_0) = f_ext
|
||||
*
|
||||
* This can be transformed to a system of 2x3n equations of first order derivative:
|
||||
*
|
||||
* x' = v
|
||||
* M v' = - C v - K (x-x_0) + f_ext
|
||||
*
|
||||
* The semi-implicit Euler scheme approximates the above with:
|
||||
*
|
||||
* x^(i+1) = x^i + \delta t * v^(i+1)
|
||||
* M v^(i+1) = M v^i + \delta t ( - C v^(i+1) - K ( x^i - x_0 ) + f^i_ext
|
||||
*
|
||||
* This is solved using implicit integration.
|
||||
*
|
||||
* @param mesh
|
||||
* @param time_step
|
||||
* @param use_stiffness_warping
|
||||
*/
|
||||
template < typename fem_mesh, typename real_type >
|
||||
inline void simulate(
|
||||
fem_mesh & mesh
|
||||
, real_type const & time_step
|
||||
, bool use_stiffness_warping,
|
||||
real_type mass_damping = 2.0,
|
||||
unsigned int min_iterations=20,
|
||||
unsigned int max_iterations=20
|
||||
)
|
||||
{
|
||||
// Some theory:
|
||||
//
|
||||
//
|
||||
// Implicit discretization of
|
||||
//
|
||||
// M d^2x/dt^2 + C dx/dt + K (x-x0) = f_ext
|
||||
//
|
||||
// Evaluate at (i+1), and use d^2x/dt^2 = (v^{i+1} - v^i)/timestep and dx/dt = v^{i+1}
|
||||
//
|
||||
// M (v^{i+1} - v^i)/timestep + C v^{i+1} + K (x^{i+1}-x0) = f_ext
|
||||
//
|
||||
// and x^{i+1} = x^i + v^{i+1}*timestep
|
||||
//
|
||||
// M (v^{i+1} - v^i)/timestep + C v^{i+1} + K ( (x^i + v^{i+1}*timestep) -x0) = f_ext
|
||||
//
|
||||
// M (v^{i+1} - v^i)/timestep + C v^{i+1} + K*x^i + timestep*K*v^{i+1} - K x0 = f_ext
|
||||
//
|
||||
// M v^{i+1} - M v^i + timestep*C v^{i+1} + timestep*timestep*K*v^{i+1} = timestep * (f_ext - K*x^i + K x0)
|
||||
//
|
||||
// (M + timestep*C + timestep*timestep*K) v^{i+1} = M v^i + timestep * (f_ext - K*x^i + K x0)
|
||||
//
|
||||
// let f0 = -K x0
|
||||
//
|
||||
// (M + timestep*C + timestep*timestep*K) v^{i+1} = M v^i + timestep * (f_ext - K*x^i -f0)
|
||||
//
|
||||
// (M + timestep*C + timestep*timestep*K) v^{i+1} = M v^i - timestep * (K*x^i + f0 - f_ext)
|
||||
//
|
||||
// so we need to solve A v^{i+1} = b for v^{i+1}, where
|
||||
//
|
||||
// A = (M + timestep*C + timestep*timestep*K)
|
||||
// b = M v^i - timestep * (K*x^i + f0 - f_ext)
|
||||
//
|
||||
// afterwards we do a position update
|
||||
//
|
||||
// x^{i+1} = x^i + v^{i+1}*timestep
|
||||
//
|
||||
// Notice that a fully implicit scheme requres K^{i+1}, however for linear elastic materials K is constant.
|
||||
|
||||
for (int n=0;n<mesh.m_nodes.size();n++)
|
||||
{
|
||||
detail::clear_stiffness_assembly_single(&mesh.m_nodes[n]);
|
||||
}
|
||||
for (int t=0;t<mesh.m_tetrahedra.size();t++)
|
||||
{
|
||||
if(use_stiffness_warping)
|
||||
detail::update_orientation_single(&mesh.m_tetrahedra[t]);
|
||||
else
|
||||
detail::reset_orientation_single(&mesh.m_tetrahedra[t]);
|
||||
}
|
||||
|
||||
for (int t=0;t<mesh.m_tetrahedra.size();t++)
|
||||
{
|
||||
detail::stiffness_assembly_single1(&mesh.m_tetrahedra[t]);
|
||||
}
|
||||
for (int i=0;i<mesh.m_tetrahedra.size();i++)
|
||||
{
|
||||
detail::add_plasticity_force_single1(mesh.m_tetrahedra[i],time_step);
|
||||
}
|
||||
|
||||
detail::dynamics_assembly(mesh,mass_damping,time_step);
|
||||
|
||||
|
||||
detail::conjugate_gradients(mesh, min_iterations, max_iterations);
|
||||
|
||||
detail::position_update(mesh,time_step);
|
||||
|
||||
}
|
||||
|
||||
} // namespace fem
|
||||
} // namespace OpenTissue
|
||||
|
||||
//OPENTISSUE_DYNAMICS_FEM_FEM_SIMULATE_H
|
||||
#endif
|
||||
@@ -0,0 +1,215 @@
|
||||
#ifndef OPENTISSUE_DYNAMICS_FEM_FEM_STIFFNESS_ASSEMBLY_H
|
||||
#define OPENTISSUE_DYNAMICS_FEM_FEM_STIFFNESS_ASSEMBLY_H
|
||||
//
|
||||
// OpenTissue Template Library
|
||||
// - A generic toolbox for physics-based modeling and simulation.
|
||||
// Copyright (C) 2008 Department of Computer Science, University of Copenhagen.
|
||||
//
|
||||
// OTTL is licensed under zlib: http://opensource.org/licenses/zlib-license.php
|
||||
//
|
||||
#include <OpenTissue/configuration.h>
|
||||
|
||||
#define p_i (&T->m_owner->m_nodes[T->m_nodes[i]])
|
||||
#define p_j (&T->m_owner->m_nodes[T->m_nodes[j]])
|
||||
|
||||
namespace OpenTissue
|
||||
{
|
||||
namespace fem
|
||||
{
|
||||
namespace detail
|
||||
{
|
||||
/**
|
||||
* 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 \in M
|
||||
* and
|
||||
* i, j \in e
|
||||
*
|
||||
* And the assembly of the force offset vector
|
||||
*
|
||||
* f0'_i = sum - Re Ke_ij x0_j
|
||||
* e \in M
|
||||
* and
|
||||
* i, j \in 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 \in M
|
||||
* and
|
||||
* i \in 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.
|
||||
*
|
||||
*
|
||||
*
|
||||
*
|
||||
*
|
||||
* @param begin
|
||||
* @param end
|
||||
*
|
||||
*/
|
||||
template<typename tetrahedron_type>
|
||||
inline void stiffness_assembly_single1(tetrahedron_type* T)
|
||||
{
|
||||
typedef typename tetrahedron_type::real_type real_type;
|
||||
typedef typename tetrahedron_type::vector3_type vector3_type;
|
||||
typedef typename tetrahedron_type::matrix3x3_type matrix3x3_type;
|
||||
|
||||
{
|
||||
matrix3x3_type & Re = T->m_Re;
|
||||
for (int i = 0; i < 4; ++i)
|
||||
{
|
||||
vector3_type f;
|
||||
f.clear();
|
||||
for (int j = 0; j < 4; ++j)
|
||||
{
|
||||
matrix3x3_type & Ke_ij = T->m_Ke[i][j];
|
||||
vector3_type & x0_j = p_j->m_model_coord;
|
||||
|
||||
f += Ke_ij * x0_j;
|
||||
if (j >= i)
|
||||
{
|
||||
|
||||
matrix3x3_type tmp = Re * Ke_ij * trans(Re);
|
||||
|
||||
p_i->K(p_j->idx()) += tmp;
|
||||
if (j > i)
|
||||
p_j->K(p_i->idx()) += trans(tmp);
|
||||
}
|
||||
}
|
||||
|
||||
p_i->m_f0 -= Re*f;
|
||||
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
} // namespace detail
|
||||
} // namespace fem
|
||||
} // namespace OpenTissue
|
||||
|
||||
//OPENTISSUE_DYNAMICS_FEM_FEM_STIFFNESS_ASSEMBLY_H
|
||||
#endif
|
||||
@@ -0,0 +1,56 @@
|
||||
#ifndef OPENTISSUE_DYNAMICS_DEFORMATION_FEM_STIFFNESSWARPING_H
|
||||
#define OPENTISSUE_DYNAMICS_DEFORMATION_FEM_STIFFNESSWARPING_H
|
||||
//
|
||||
// OpenTissue Template Library
|
||||
// - A generic toolbox for physics-based modeling and simulation.
|
||||
// Copyright (C) 2008 Department of Computer Science, University of Copenhagen.
|
||||
//
|
||||
// OTTL is licensed under zlib: http://opensource.org/licenses/zlib-license.php
|
||||
//
|
||||
#include <OpenTissue/configuration.h>
|
||||
|
||||
namespace OpenTissue
|
||||
{
|
||||
namespace fem
|
||||
{
|
||||
namespace detail
|
||||
{
|
||||
template <typename math_types>
|
||||
class TetrahedronTraits
|
||||
{
|
||||
public:
|
||||
|
||||
typedef typename math_types::real_type real_type;
|
||||
typedef typename math_types::vector3_type vector3_type;
|
||||
typedef typename math_types::matrix3x3_type matrix3x3_type;
|
||||
|
||||
public:
|
||||
|
||||
real_type m_young;
|
||||
real_type m_poisson;
|
||||
real_type m_density;
|
||||
|
||||
matrix3x3_type m_Ke[4][4]; ///< Stiffness element matrix
|
||||
matrix3x3_type m_Re; ///< Rotational warp of tetrahedron.
|
||||
real_type m_V; ///< Volume of tetrahedron
|
||||
|
||||
vector3_type m_e10; ///< edge from p0 to p1
|
||||
vector3_type m_e20; ///< edge from p0 to p2
|
||||
vector3_type m_e30; ///< edge from p0 to p3
|
||||
|
||||
//--- Stuff used exclusive by plastic effects
|
||||
|
||||
vector3_type m_B[4]; ///< placeholders for Jacobian of shapefunctions: B = SN.
|
||||
vector3_type m_D; ///< Elasticity Matrix in vector from
|
||||
real_type m_plastic[6]; ///< Plastic strain tensor.
|
||||
real_type m_yield;
|
||||
real_type m_creep;
|
||||
real_type m_max;
|
||||
};
|
||||
|
||||
} // namespace detail
|
||||
} // namespace fem
|
||||
} // namespace OpenTissue
|
||||
|
||||
//OPENTISSUE_DYNAMICS_DEFORMATION_FEM_STIFFNESSWARPING_H
|
||||
#endif
|
||||
@@ -0,0 +1,38 @@
|
||||
#ifndef OPENTISSUE_DYNAMICS_FEM_FEM_UNIFORM_DENSITY_H
|
||||
#define OPENTISSUE_DYNAMICS_FEM_FEM_UNIFORM_DENSITY_H
|
||||
//
|
||||
// OpenTissue Template Library
|
||||
// - A generic toolbox for physics-based modeling and simulation.
|
||||
// Copyright (C) 2008 Department of Computer Science, University of Copenhagen.
|
||||
//
|
||||
// OTTL is licensed under zlib: http://opensource.org/licenses/zlib-license.php
|
||||
//
|
||||
#include <OpenTissue/configuration.h>
|
||||
|
||||
namespace OpenTissue
|
||||
{
|
||||
namespace fem
|
||||
{
|
||||
namespace detail
|
||||
{
|
||||
/**
|
||||
* Set Uniform Density.
|
||||
*
|
||||
* @param begin
|
||||
* @param end
|
||||
* @param density
|
||||
*/
|
||||
template<typename real_type, typename tetrahedron_iterator>
|
||||
inline void uniform_density(tetrahedron_iterator const & begin, tetrahedron_iterator const & end,real_type const & density)
|
||||
{
|
||||
assert(density>0 || !"uniform_density(): density must be positive");
|
||||
for (tetrahedron_iterator T = begin;T!=end;++T)
|
||||
T->m_density = density;
|
||||
}
|
||||
|
||||
} // namespace detail
|
||||
} // namespace fem
|
||||
} // namespace OpenTissue
|
||||
|
||||
//OPENTISSUE_DYNAMICS_FEM_FEM_UNIFORM_DENSITY_H
|
||||
#endif
|
||||
@@ -0,0 +1,44 @@
|
||||
#ifndef OPENTISSUE_DYNAMICS_FEM_FEM_UNIFORM_POISSON_H
|
||||
#define OPENTISSUE_DYNAMICS_FEM_FEM_UNIFORM_POISSON_H
|
||||
//
|
||||
// OpenTissue Template Library
|
||||
// - A generic toolbox for physics-based modeling and simulation.
|
||||
// Copyright (C) 2008 Department of Computer Science, University of Copenhagen.
|
||||
//
|
||||
// OTTL is licensed under zlib: http://opensource.org/licenses/zlib-license.php
|
||||
//
|
||||
#include <OpenTissue/configuration.h>
|
||||
|
||||
namespace OpenTissue
|
||||
{
|
||||
namespace fem
|
||||
{
|
||||
namespace detail
|
||||
{
|
||||
|
||||
/**
|
||||
*
|
||||
* @param begin
|
||||
* @param end
|
||||
* @param poisson
|
||||
*/
|
||||
template<typename real_type, typename tetrahedron_iterator>
|
||||
inline void uniform_poisson(
|
||||
tetrahedron_iterator const & begin
|
||||
, tetrahedron_iterator const & end
|
||||
, real_type const & poisson
|
||||
)
|
||||
{
|
||||
assert(poisson>=0 || !"uniform_poisson() : Poisson ratio must be non-negative");
|
||||
assert(poisson<=0.5 || !"uniform_poisson() : Poisson ratio must not be larger than a half");
|
||||
|
||||
for (tetrahedron_iterator T = begin;T!=end;++T)
|
||||
T->m_poisson = poisson;
|
||||
}
|
||||
|
||||
} // namespace detail
|
||||
} // namespace fem
|
||||
} // namespace OpenTissue
|
||||
|
||||
//OPENTISSUE_DYNAMICS_FEM_FEM_UNIFORM_POISSON_H
|
||||
#endif
|
||||
@@ -0,0 +1,43 @@
|
||||
#ifndef OPENTISSUE_DYNAMICS_FEM_FEM_UNIFORM_YOUNG_H
|
||||
#define OPENTISSUE_DYNAMICS_FEM_FEM_UNIFORM_YOUNG_H
|
||||
//
|
||||
// OpenTissue Template Library
|
||||
// - A generic toolbox for physics-based modeling and simulation.
|
||||
// Copyright (C) 2008 Department of Computer Science, University of Copenhagen.
|
||||
//
|
||||
// OTTL is licensed under zlib: http://opensource.org/licenses/zlib-license.php
|
||||
//
|
||||
#include <OpenTissue/configuration.h>
|
||||
|
||||
namespace OpenTissue
|
||||
{
|
||||
namespace fem
|
||||
{
|
||||
namespace detail
|
||||
{
|
||||
|
||||
/**
|
||||
*
|
||||
* @param begin
|
||||
* @param end
|
||||
* @param young
|
||||
*/
|
||||
template<typename real_type, typename tetrahedron_iterator>
|
||||
inline void uniform_young(
|
||||
tetrahedron_iterator const & begin
|
||||
, tetrahedron_iterator const & end
|
||||
,real_type const & young
|
||||
)
|
||||
{
|
||||
assert(young>=0 || !"uniform_young(): Young modulus must not be negative");
|
||||
|
||||
for (tetrahedron_iterator T = begin;T!=end;++T)
|
||||
T->m_young = young;
|
||||
}
|
||||
|
||||
} // namespace detail
|
||||
} // namespace fem
|
||||
} // namespace OpenTissue
|
||||
|
||||
//OPENTISSUE_DYNAMICS_FEM_FEM_UNIFORM_YOUNG_H
|
||||
#endif
|
||||
@@ -0,0 +1,222 @@
|
||||
#ifndef OPENTISSUE_DYNAMICS_FEM_FEM_UPDATE_ORIENTATION_H
|
||||
#define OPENTISSUE_DYNAMICS_FEM_FEM_UPDATE_ORIENTATION_H
|
||||
//
|
||||
// OpenTissue Template Library
|
||||
// - A generic toolbox for physics-based modeling and simulation.
|
||||
// Copyright (C) 2008 Department of Computer Science, University of Copenhagen.
|
||||
//
|
||||
// OTTL is licensed under zlib: http://opensource.org/licenses/zlib-license.php
|
||||
//
|
||||
#include <OpenTissue/configuration.h>
|
||||
|
||||
namespace OpenTissue
|
||||
{
|
||||
namespace fem
|
||||
{
|
||||
namespace detail
|
||||
{
|
||||
|
||||
/**
|
||||
*
|
||||
*/
|
||||
template<typename tetrahedron_type>
|
||||
inline void update_orientation_single(tetrahedron_type* T)
|
||||
{
|
||||
typedef typename tetrahedron_type::real_type real_type;
|
||||
typedef typename tetrahedron_type::vector3_type vector3_type;
|
||||
typedef typename tetrahedron_type::matrix3x3_type matrix3x3_type;
|
||||
|
||||
{
|
||||
real_type div6V = 1.0 / T->m_V*6.0;
|
||||
//--- The derivation in the orignal paper on stiffness warping were stated as
|
||||
//---
|
||||
//--- | n1^T |
|
||||
//--- N = | n2^T | : WCS -> U
|
||||
//--- | n3^T |
|
||||
//---
|
||||
//--- | n1'^T |
|
||||
//--- N' = | n2'^T | : WCS -> D
|
||||
//--- | n3'^T |
|
||||
//---
|
||||
//--- From which we have
|
||||
//---
|
||||
//--- R = N' N^T
|
||||
//---
|
||||
//--- This is valid under the assumption that the n-vectors form a orthonormal basis. In that
|
||||
//--- paper the n-vectors were determined using a heuristic approach. Besides
|
||||
//--- the rotation were computed on a per vertex basis not on a per-tetrahedron
|
||||
//--- bases as we have outline above.
|
||||
//---
|
||||
//---
|
||||
//--- Later Muller et. al. used barycentric coordinates to find the transform. We
|
||||
//--- will here go into details on this method. Let the deformed corners be q0, q1,
|
||||
//--- q2, and q3 and the undeformed corners p0, p1, p2,and p3. Looking at some
|
||||
//--- point p = [x y z 1]^T inside the undeformed tetrahedron this can be written
|
||||
//---
|
||||
//--- | w0 |
|
||||
//--- p = | p0 p1 p2 p3 | | w1 | = P w (*1)
|
||||
//--- | 1 1 1 1 | | w2 |
|
||||
//--- | w3 |
|
||||
//---
|
||||
//--- The same point in the deformed tetrahedron has the same barycentric coordinates
|
||||
//--- which mean
|
||||
//---
|
||||
//--- | w0 |
|
||||
//--- q = | q0 q1 q2 q3 | | w1 | = Q w (*2)
|
||||
//--- | 1 1 1 1 | | w2 |
|
||||
//--- | w3 |
|
||||
//---
|
||||
//--- We can now use (*1) to solve for w and insert this into (*2), this yields
|
||||
//---
|
||||
//--- q = Q P^{-1} p
|
||||
//---
|
||||
//--- The matirx Q P^{-1} transforms p into q. Due to P and Q having their fourth rows
|
||||
//--- equal to 1 it can be shown that this matrix have the block structure
|
||||
//---
|
||||
//--- Q P^{-1} = | R t | (*3)
|
||||
//--- | 0^T 1 |
|
||||
//---
|
||||
//--- Which we recognize as a transformation matrix of homegeneous coordinates. The t vector
|
||||
//--- gives the translation, the R matrix includes, scaling, shearing and rotation.
|
||||
//---
|
||||
//--- We thus need to extract the rotational information from this R matrix. In their
|
||||
//--- paper Mueller et. al. suggest using Polar Decompostions (cite Shoemake and Duff; Etzmuss).
|
||||
//--- However, a simpler although more imprecise approach would simply be to apply a Grahram-Schimdt
|
||||
//--- orthonormalization to transform R into an orthonormal matrix. This seems to work quite well
|
||||
//--- in practice. We have not observed any visual difference in using polar decomposition or
|
||||
//--- orthonormalization.
|
||||
//---
|
||||
//--- Now for some optimizations. Since we are only interested in the R part of (*3) we can
|
||||
//--- compute this more efficiently exploiting the fact that barycentric coordinates sum
|
||||
//--- up to one. If we substitute w0 = 1 - w1 - w2 - w3 in (*1) and (*2) we can throw away
|
||||
//--- the fourth rows since they simply state 1=1, which is trivially true.
|
||||
//---
|
||||
//--- | 1-w1-w2-w3 |
|
||||
//--- p = | p0 p1 p2 p3 | | w1 |
|
||||
//--- | w2 |
|
||||
//--- | w3 |
|
||||
//---
|
||||
//--- Which is
|
||||
//---
|
||||
//--- | 1 |
|
||||
//--- p = | p0 (p1-p0) (p2-p0) (p3-p0) | | w1 |
|
||||
//--- | w2 |
|
||||
//--- | w3 |
|
||||
//---
|
||||
//--- We can move the first column over on the left hand side
|
||||
//---
|
||||
//--- | w1 |
|
||||
//--- (p-p0) = | (p1-p0) (p2-p0) (p3-p0) | | w2 |
|
||||
//--- | w3 |
|
||||
//---
|
||||
//--- Introducing e10 = p1-p0, e20 = p2-p0, e30 = p3-p0, and E = [e10 e20 e30 ] we have
|
||||
//---
|
||||
//--- w1
|
||||
//--- (p-p0) = E w2 (*5)
|
||||
//--- w3
|
||||
//---
|
||||
//--- Similar for *(2) we have
|
||||
//---
|
||||
//--- w1
|
||||
//--- (q-q0) = E' w2 (*6)
|
||||
//--- w3
|
||||
//---
|
||||
//--- Where E' = [e10' e20' e3'] and e10' = q1-q0, e20' = q2-q0, e30' = q3-q0. Now
|
||||
//--- inverting (*5) and insertion into (*6) yields
|
||||
//---
|
||||
//--- (q-q0) = E' E^{-1} (p-p0)
|
||||
//---
|
||||
//--- By comparison with (*3) we see that
|
||||
//---
|
||||
//--- R = E' E^{-1}
|
||||
//---
|
||||
//--- Using Cramers rule the inverse of E can be written as
|
||||
//---
|
||||
//--- | (e2 x e3)^T |
|
||||
//--- E^{-1} = 1/6V | (e3 x e1)^T |
|
||||
//--- | (e1 x e2)^T |
|
||||
//---
|
||||
//--- This can easily be confirmed by straigthforward computation
|
||||
//---
|
||||
//--- | e1 \cdot (e2 x e3) e2 \cdot (e2 x e3) e3 \cdot (e2 x e3) |
|
||||
//--- E^{-1} E = 1/6v | e1 \cdot (e3 x e1) e2 \cdot (e3 x e1) e3 \cdot (e3 x e1) |
|
||||
//--- | e1 \cdot (e1 x e2) e2 \cdot (e1 x e2) e3 \cdot (e1 x e2) |
|
||||
//---
|
||||
//--- | 6V 0 0 |
|
||||
//--- = 1/6V | 0 6V 0 | = I
|
||||
//--- | 0 0 6V |
|
||||
//---
|
||||
//--- Using the notation
|
||||
//---
|
||||
//--- n1 = e2 x e3 / 6V
|
||||
//--- n2 = e3 x e1 / 6V
|
||||
//--- n3 = e1 x e2 / 6V
|
||||
//---
|
||||
//--- we write
|
||||
//---
|
||||
//--- | n1^T |
|
||||
//--- E^{-1} = | n2^T |
|
||||
//--- | n3^T |
|
||||
//---
|
||||
//--- And we end up with
|
||||
//---
|
||||
//--- | n1^T |
|
||||
//--- R = [e10' e20' e30'] | n2^T |
|
||||
//--- | n3^T |
|
||||
//---
|
||||
//--- Observe that E' is very in-expensive to compute and all non-primed quantities (n1, n2, n3) can
|
||||
//--- be precomputed and stored on a per tetrahedron basis if there is enough memory available. Even
|
||||
//--- in case where memory is not available n1, n2 and n3 are quite cheap to compute.
|
||||
//---
|
||||
|
||||
real_type e1x = T->m_e10(0);
|
||||
real_type e1y = T->m_e10(1);
|
||||
real_type e1z = T->m_e10(2);
|
||||
real_type e2x = T->m_e20(0);
|
||||
real_type e2y = T->m_e20(1);
|
||||
real_type e2z = T->m_e20(2);
|
||||
real_type e3x = T->m_e30(0);
|
||||
real_type e3y = T->m_e30(1);
|
||||
real_type e3z = T->m_e30(2);
|
||||
real_type n1x = (e2y * e3z - e3y * e2z) * div6V;
|
||||
real_type n1y = (e3x * e2z - e2x * e3z) * div6V;
|
||||
real_type n1z = (e2x * e3y - e3x * e2y) * div6V;
|
||||
real_type n2x = (e1z * e3y - e1y * e3z) * div6V;
|
||||
real_type n2y = (e1x * e3z - e1z * e3x) * div6V;
|
||||
real_type n2z = (e1y * e3x - e1x * e3y) * div6V;
|
||||
real_type n3x = (e1y * e2z - e1z * e2y) * div6V;
|
||||
real_type n3y = (e1z * e2x - e1x * e2z) * div6V;
|
||||
real_type n3z = (e1x * e2y - e1y * e2x) * div6V;
|
||||
const vector3_type & p0 = T->m_owner->m_nodes[T->m_nodes[0]].m_coord;
|
||||
const vector3_type & p1 = T->m_owner->m_nodes[T->m_nodes[1]].m_coord;
|
||||
const vector3_type & p2 = T->m_owner->m_nodes[T->m_nodes[2]].m_coord;
|
||||
const vector3_type & p3 = T->m_owner->m_nodes[T->m_nodes[3]].m_coord;
|
||||
|
||||
e1x = p1(0) - p0(0);
|
||||
e1y = p1(1) - p0(1);
|
||||
e1z = p1(2) - p0(2);
|
||||
e2x = p2(0) - p0(0);
|
||||
e2y = p2(1) - p0(1);
|
||||
e2z = p2(2) - p0(2);
|
||||
e3x = p3(0) - p0(0);
|
||||
e3y = p3(1) - p0(1);
|
||||
e3z = p3(2) - p0(2);
|
||||
T->m_Re(0,0) = e1x * n1x + e2x * n2x + e3x * n3x; T->m_Re(0,1) = e1x * n1y + e2x * n2y + e3x * n3y; T->m_Re(0,2) = e1x * n1z + e2x * n2z + e3x * n3z;
|
||||
T->m_Re(1,0) = e1y * n1x + e2y * n2x + e3y * n3x; T->m_Re(1,1) = e1y * n1y + e2y * n2y + e3y * n3y; T->m_Re(1,2) = e1y * n1z + e2y * n2z + e3y * n3z;
|
||||
T->m_Re(2,0) = e1z * n1x + e2z * n2x + e3z * n3x; T->m_Re(2,1) = e1z * n1y + e2z * n2y + e3z * n3y; T->m_Re(2,2) = e1z * n1z + e2z * n2z + e3z * n3z;
|
||||
|
||||
T->m_Re = ortonormalize(T->m_Re);
|
||||
|
||||
//matrix3x3_type M = T->m_Re,S;
|
||||
//OpenTissue::PolarDecomposition3x3 decomp;
|
||||
//decomp.eigen(M, T->m_Re, S);//--- Etzmuss-style
|
||||
//decomp.newton(M, T->m_Re );//--- Shoemake-Duff-style
|
||||
}
|
||||
}
|
||||
|
||||
} // namespace detail
|
||||
} // namespace fem
|
||||
} // namespace OpenTissue
|
||||
|
||||
//OPENTISSUE_DYNAMICS_FEM_FEM_UPDATE_ORIENTATION_H
|
||||
#endif
|
||||
Reference in New Issue
Block a user