00001 #ifndef OPENTISSUE_DYNAMICS_FEM_FEM_COMPUTE_MASS_H 00002 #define OPENTISSUE_DYNAMICS_FEM_FEM_COMPUTE_MASS_H 00003 // 00004 // OpenTissue Template Library 00005 // - A generic toolbox for physics-based modeling and simulation. 00006 // Copyright (C) 2008 Department of Computer Science, University of Copenhagen. 00007 // 00008 // OTTL is licensed under zlib: http://opensource.org/licenses/zlib-license.php 00009 // 00010 #include <OpenTissue/configuration.h> 00011 00012 #include <limits> 00013 00014 namespace OpenTissue 00015 { 00016 namespace fem 00017 { 00018 namespace detail 00019 { 00020 00030 template<typename fem_mesh> 00031 inline void compute_mass(fem_mesh & mesh) 00032 { 00033 typedef typename fem_mesh::node_iterator node_iterator; 00034 typedef typename fem_mesh::tetrahedron_iterator tetrahedron_iterator; 00035 typedef typename fem_mesh::real_type real_type; 00036 00037 //--- 00038 //--- A mass matrix is a discrete representation of a continuous mass distribution. 00039 //--- To compute our mass matrix for a tetrahedral element with linear shape 00040 //--- functions we need the formula (pp. 266 in Cook) 00041 //--- 00042 //--- a!b!c!d! 00043 //--- \int_V N_1^a N_2^b N_3^c N_4^d dV = 6V -------------------------- (**) 00044 //--- (3 + a + b +c + d)! 00045 //--- 00046 //--- A consistent element mass matrix (pp. 376 Cook) is defined as 00047 //--- 00048 //--- m = \int_V \rho N^T N dV (***) 00049 //--- 00050 //--- This equation can be derived from work balance, the details of which is unimportant 00051 //--- here (look Cook pp. 375-376 for details). 00052 //--- Assumping \rho is constant over each tetrahedral element and using the linear shape 00053 //--- functions the above definition (***) results in 00054 //--- 00055 //--- |N_1| 00056 //--- m = \rho \int_V |N_2| |N_1 N_2 N_3 N_4| dV 00057 //--- |N_3| 00058 //--- |N_4| 00059 //--- 00060 //--- |(N_1 N_1) (N_1 N_2) (N_1 N_3) (N_1 N_4)| 00061 //--- m = \rho \int_V |(N_2 N_1) (N_2 N_2) (N_2 N_3) (N_2 N_4)| dV 00062 //--- |(N_3 N_1) (N_3 N_2) (N_3 N_3) (N_3 N_4)| 00063 //--- |(N_4 N_1) (N_4 N_2) (N_4 N_3) (N_4 N_4)| 00064 //--- 00065 //--- by (**) 00066 //--- 00067 //--- | 2 1 1 1| 00068 //--- m = \rho V/20 | 1 2 1 1| (****) 00069 //--- | 1 1 2 1| 00070 //--- | 1 1 1 2| 00071 //--- 00072 //--- V 00073 //--- m_ij = \rho --- (1+delta_ij) 00074 //--- 20 00075 //--- 00076 //--- in 3D this means that for the tetrahedral element 00077 //--- 00078 //--- | 2 2 2 1 1 1 1 1 1 1 1 1 | 00079 //--- | 2 2 2 1 1 1 1 1 1 1 1 1 | 00080 //--- | 2 2 2 1 1 1 1 1 1 1 1 1 | 00081 //--- | | 00082 //--- | 1 1 1 2 2 2 1 1 1 1 1 1 | 00083 //--- | 1 1 1 2 2 2 1 1 1 1 1 1 | 00084 //--- V | 1 1 1 2 2 2 1 1 1 1 1 1 | 00085 //--- Me = \rho --- | | 00086 //--- 20 | 1 1 1 1 1 1 2 2 2 1 1 1 | 00087 //--- | 1 1 1 1 1 1 2 2 2 1 1 1 | 00088 //--- | 1 1 1 1 1 1 2 2 2 1 1 1 | 00089 //--- | | 00090 //--- | 1 1 1 1 1 1 1 1 1 2 2 2 | 00091 //--- | 1 1 1 1 1 1 1 1 1 2 2 2 | 00092 //--- | 1 1 1 1 1 1 1 1 1 2 2 2 | 00093 //--- 00094 //--- Notice that in order to obtain the global/system mass matrix an assembly similar to the 00095 //--- stiffnees matrix assembly must be carried out. Further, the global M matrix will 00096 //--- have the same sub-block pattern as the global K matrix. 00097 //--- 00098 //--- A consistent mass matrix is often not used in computer graphics. Instead and 00099 //--- ad-hoc approach named ``lumped'' mass matrix is applied. 00100 //--- The lumped mass matrix is obtained by placing particle masses at the nodes. 00101 //--- This corresponds to shifting all the masses in the rows of (****) onto the 00102 //--- diagonal. In 3D this yields the element mass matrix 00103 //--- 00104 //--- | 1 0 0 0 0 0 0 0 0 0 0 0 | 00105 //--- | 0 1 0 0 0 0 0 0 0 0 0 0 | 00106 //--- | 0 0 1 0 0 0 0 0 0 0 0 0 | 00107 //--- | | 00108 //--- | 0 0 0 1 0 0 0 0 0 0 0 0 | 00109 //--- | 0 0 0 0 1 0 0 0 0 0 0 0 | 00110 //--- V | 0 0 0 0 0 1 0 0 0 0 0 0 | 00111 //--- Me = \rho --- | | 00112 //--- 4 | 0 0 0 0 0 0 1 0 0 0 0 0 | 00113 //--- | 0 0 0 0 0 0 0 1 0 0 0 0 | 00114 //--- | 0 0 0 0 0 0 0 0 1 0 0 0 | 00115 //--- | | 00116 //--- | 0 0 0 0 0 0 0 0 0 1 0 0 | 00117 //--- | 0 0 0 0 0 0 0 0 0 0 1 0 | 00118 //--- | 0 0 0 0 0 0 0 0 0 0 0 1 | 00119 //--- 00120 //--- Thus a lumped mass matrix is diagonal whereas a consistent mass matrix 00121 //--- is not. Observe that the global mass matrix would also diagonal and the 00122 //--- assembly is simplified to an iteration over all tetrahedra, while 00123 //--- incementing the nodal mass by one fourth of the tetrahedral mass. 00124 //--- 00125 //--- for each node n 00126 //--- mass(n) = 0 00127 //--- next n 00128 //--- for each tetrahedron e 00129 //--- mass(n_i) += \rho_e Ve / 4 00130 //--- mass(n_j) += \rho_e Ve / 4 00131 //--- mass(n_k) += \rho_e Ve / 4 00132 //--- mass(n_m) += \rho_e Ve / 4 00133 //--- next e 00134 //--- 00135 //--- where n_i,n_j,n_k and n_m are the four nodes of the e'th tetrahedron. 00136 //--- 00137 //--- The advantage of lumping is less storage and higher performace. On the downside 00138 //--- lumping introduces a discontinouty in the displacement field. 00139 //--- 00140 //--- Obrien.shen state that the errors in lumping is negligeble for small-size course 00141 //--- meshes used in computer graphics. However, for finer meshes the errors becomes 00142 //--- noticeable. 00143 //--- 00144 //--- There do exist other approaches for computing mass matrices, even mehtods which 00145 //--- combine other methods. We refer the interested reader to Cook for more details. Here 00146 //--- we have limited our selfes to the two most common methods. 00147 //--- 00148 //--- It is worthwhile to notice that under the reasonable assumptions that V and \rho are 00149 //--- positive for all elements both the element mass matrices and the global mass matrices 00150 //--- are symmetric positive definite matrices. 00151 00152 node_iterator Nbegin = mesh.node_begin(); 00153 node_iterator Nend = mesh.node_end(); 00154 for (node_iterator N = Nbegin;N!=Nend;++N) 00155 { 00156 if(N->m_fixed) 00157 N->m_mass = std::numeric_limits<real_type>::max(); 00158 else 00159 N->m_mass = 0; 00160 } 00161 00162 tetrahedron_iterator Tbegin = mesh.tetrahedron_begin(); 00163 tetrahedron_iterator Tend = mesh.tetrahedron_end(); 00164 for (tetrahedron_iterator T = Tbegin;T!=Tend;++T) 00165 { 00166 real_type amount = T->m_density * T->m_V *.25; 00167 T->i()->m_mass += amount; 00168 T->j()->m_mass += amount; 00169 T->k()->m_mass += amount; 00170 T->m()->m_mass += amount; 00171 } 00172 } 00173 00174 } // namespace detail 00175 } // namespace fem 00176 } // namespace OpenTissue 00177 00178 //OPENTISSUE_DYNAMICS_FEM_FEM_COMPUTE_MASS_H 00179 #endif
1.7.1