00001 #ifndef OPENTISSUE_DYNAMICS_VERSATILE_VERSATILE_AREA_CONSTRAINT_H 00002 #define OPENTISSUE_DYNAMICS_VERSATILE_VERSATILE_AREA_CONSTRAINT_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 namespace OpenTissue 00013 { 00014 namespace versatile 00015 { 00016 namespace detail 00017 { 00018 template <typename versatile_types> 00019 class AreaConstraint 00020 { 00021 public: 00022 00023 typedef typename versatile_types::value_traits value_traits; 00024 typedef typename versatile_types::real_type real_type; 00025 typedef typename versatile_types::vector3_type vector3_type; 00026 typedef typename versatile_types::node_type node_type; 00027 00028 public: 00029 00030 real_type m_k; 00031 real_type m_b; 00032 real_type m_A0; 00033 node_type * m_ni; 00034 node_type * m_nj; 00035 node_type * m_nk; 00036 00037 public: 00038 00039 void initialize(node_type & ni, node_type & nj, node_type & nk) 00040 { 00041 using std::sqrt; 00042 00043 m_ni = & ni; 00044 m_nj = & nj; 00045 m_nk = & nk; 00046 vector3_type A = cross( (m_nj->m_coord - m_ni->m_coord) , (m_nk->m_coord - m_ni->m_coord) ); 00047 m_A0 = value_traits::half()*sqrt(A*A); 00048 } 00049 00050 void apply() const 00051 { 00052 using std::sqrt; 00053 00054 vector3_type e_jk = (m_nj->m_coord - m_nk->m_coord); 00055 vector3_type e_ki = (m_nk->m_coord - m_ni->m_coord); 00056 vector3_type e_ij = (m_ni->m_coord - m_nj->m_coord); 00057 00058 vector3_type A = cross(e_ki , e_ij); 00059 real_type A_norm = sqrt(A*A); 00060 real_type C = ( (value_traits::half()*A_norm - m_A0) / m_A0 ); 00061 real_type nominator = value_traits::one() / (value_traits::two()*m_A0 * A_norm); 00062 vector3_type dCdi = cross( e_jk , A )*nominator; 00063 vector3_type dCdj = cross( e_ki , A )*nominator; 00064 vector3_type dCdk = cross( e_ij , A )*nominator; 00065 00066 real_type factor = value_traits::zero(); 00067 00068 if(m_k>value_traits::zero()) 00069 factor -= m_k *C; 00070 00071 if(m_b>value_traits::zero()) 00072 factor -= m_b * (dCdi*m_ni->m_v + dCdj*m_nj->m_v + dCdk*m_nk->m_v); 00073 00074 if(factor) 00075 { 00076 m_ni->m_f_con += factor*dCdi; 00077 m_nj->m_f_con += factor*dCdj; 00078 m_nk->m_f_con += factor*dCdk; 00079 } 00080 } 00081 00082 real_type compute_internal_energy() 00083 { 00084 vector3_type e_ki = (m_nk->m_coord - m_ni->m_coord); 00085 vector3_type e_ij = (m_ni->m_coord - m_nj->m_coord); 00086 vector3_type A = cross(e_ki , e_ij); 00087 real_type A_norm = sqrt(A*A); 00088 real_type C = ( (value_traits::half()*A_norm-m_A0) / m_A0 ); 00089 return value_traits::half()*m_k*C*C; 00090 } 00091 }; 00092 00093 } // namespace detail 00094 } // namespace versatile 00095 } // namespace OpenTissue 00096 00097 //OPENTISSUE_DYNAMICS_VERSATILE_VERSATILE_AREA_CONSTRAINT_H 00098 #endif
1.7.1