HUMIM Tracker: /home/hauberg/Dokumenter/Capture/humim-tracker-0.1/src/OpenTissue/OpenTissue/dynamics/fem/fem_compute

Go to the documentation of this file.
00001 #ifndef OPENTISSUE_DYNAMICS_FEM_FEM_COMPUTE_KE_H
00002 #define OPENTISSUE_DYNAMICS_FEM_FEM_COMPUTE_KE_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 fem
00015   {
00016     namespace detail
00017     {
00029       template<typename real_type, typename vector3_type, typename matrix3x3_type>
00030       inline void compute_Ke(
00031         vector3_type const & p0,
00032         vector3_type const & p1,
00033         vector3_type const & p2,
00034         vector3_type const & p3,
00035         real_type const & E,
00036         real_type const & nu,
00037         matrix3x3_type Ke[4][4]
00038       )
00039       {
00040         using std::fabs;
00041 
00042         real_type d = p0(0);
00043         real_type d4 = p0(1);
00044         real_type d8 = p0(2);
00045         real_type d1 = p1(0) - d;
00046         real_type d5 = p1(1) - d4;
00047         real_type d9 = p1(2) - d8;
00048         real_type d2 = p2(0) - d;
00049         real_type d6 = p2(1) - d4;
00050         real_type d10 = p2(2) - d8;
00051         real_type d3 = p3(0) - d;
00052         real_type d7 = p3(1) - d4;
00053         real_type d11 = p3(2) - d8;
00054         real_type d12 = (d1 * d6 * d11 + d2 * d7 * d9 + d3 * d5 * d10) - d1 * d7 * d10 - d2 * d5 * d11 - d3 * d6 * d9;
00055         real_type d13 = 1.0 / d12;
00056         real_type d14 = fabs(d12 / 6.0);
00057         vector3_type B[4];
00058         B[0].clear();
00059         B[1].clear();
00060         B[2].clear();
00061         B[3].clear();
00062         B[1](0) = (d10 * d7 - d6 * d11) * d13;
00063         B[2](0) = (d5 * d11 - d9 * d7) * d13;
00064         B[3](0) = (d9 * d6 - d5 * d10) * d13;
00065         B[0](0) = -B[1](0) - B[2](0) - B[3](0);
00066         B[1](1) = (d2 * d11 - d10 * d3) * d13;
00067         B[2](1) = (d9 * d3 - d1 * d11) * d13;
00068         B[3](1) = (d1 * d10 - d9 * d2) * d13;
00069         B[0](1) = -B[1](1) - B[2](1) - B[3](1);
00070         B[1](2) = (d6 * d3 - d2 * d7) * d13;
00071         B[2](2) = (d1 * d7 - d5 * d3) * d13;
00072         B[3](2) = (d5 * d2 - d1 * d6) * d13;
00073         B[0](2) = -B[1](2) - B[2](2) - B[3](2);
00074         real_type d15 = E / (1.0 + nu) / (1.0 - 2 * nu);
00075         real_type d16 = (1.0 - nu) * d15;
00076         real_type d17 = nu * d15;
00077         real_type d18 = E / 2 / (1.0 + nu);
00078         for(int i = 0; i < 4; ++i)
00079         {
00080           for(int j = 0; j < 4; ++j)
00081           {
00082             real_type d19 = B[i](0);
00083             real_type d20 = B[i](1);
00084             real_type d21 = B[i](2);
00085             real_type d22 = B[j](0);
00086             real_type d23 = B[j](1);
00087             real_type d24 = B[j](2);
00088             Ke[i][j](0,0) = d16 * d19 * d22 + d18 * (d20 * d23 + d21 * d24);
00089             Ke[i][j](0,1) = d17 * d19 * d23 + d18 * (d20 * d22);
00090             Ke[i][j](0,2) = d17 * d19 * d24 + d18 * (d21 * d22);
00091             Ke[i][j](1,0) = d17 * d20 * d22 + d18 * (d19 * d23);
00092             Ke[i][j](1,1) = d16 * d20 * d23 + d18 * (d19 * d22 + d21 * d24);
00093             Ke[i][j](1,2) = d17 * d20 * d24 + d18 * (d21 * d23);
00094             Ke[i][j](2,0) = d17 * d21 * d22 + d18 * (d19 * d24);
00095             Ke[i][j](2,1) = d17 * d21 * d23 + d18 * (d20 * d24);
00096             Ke[i][j](2,2) = d16 * d21 * d24 + d18 * (d20 * d23 + d19 * d22);
00097             Ke[i][j] *= d14;
00098           }
00099         }
00100       }
00101 
00106       template<typename real_type, typename vector3_type, typename matrix3x3_type>
00107       inline void compute_Ke(
00108         vector3_type * B,
00109         vector3_type const& D,
00110         real_type const & volume,
00111         matrix3x3_type Ke[4][4]
00112       )
00113       {
00114         for (unsigned int i=0; i<4; ++i)
00115           for (unsigned int j=i; j<4; ++j)
00116             compute_Ke_ij( B[i], D, B[j], volume, Ke[i][j] );
00117 
00118         for (unsigned int i=1; i<4; ++i)
00119           for (unsigned int j=0; j<i; ++j)
00120             Ke[i][j] = trans(Ke[j][i]);
00121       }
00122 
00132       template<typename real_type, typename vector3_type, typename matrix3x3_type>
00133       inline void compute_Ke_ij(
00134         vector3_type const& Bi,
00135         vector3_type const& D,
00136         vector3_type const& Bj,
00137         real_type const & volume,
00138         matrix3x3_type & Ke_ij
00139         )
00140       {
00141         assert(Ke_ij.size1() == 3 || !"compute_Ke_ij(): Incompatible dimension");
00142         assert(Ke_ij.size2() == 3 || !"compute_Ke_ij(): Incompatible dimension");
00143 
00144         //
00145         // Calculate Ke_ij = Bi^T D Bj V
00146         //
00147         // The matrix Bi is given by:
00148         //
00149         //          | bi 0  0  |
00150         //          | 0  ci 0  |
00151         //   Bi =   | 0  0  di |
00152         //          | ci bi 0  |
00153         //          | di 0  bi |
00154         //          | 0  di ci |
00155         //
00156         //  Similar for Bj. The elasticity matrix D is given by
00157         //
00158         //
00159         //                  -                                             -     -                        -
00160         //                  | 1-nu  nu   nu      0         0        0     |     | D0  D1  D1  0   0   0  |
00161         //       young      |  nu  1-nu  nu      0         0        0     |     | D1  D0  D1  0   0   0  |
00162         // D= ------------- |  nu   nu  1-nu     0         0        0     |     | D1  D1  D0  0   0   0  |
00163         //    (1+nu)(1-2nu) |  0    0   0    (1-2nu)/2     0        0     |  =  | 0   0   0   D2  0   0  |
00164         //                  |  0    0   0        0     (1-2nu)/2    0     |     | 0   0   0   0   D2  0  |
00165         //                  |  0    0   0        0         0    (1-2nu)/2 |     | 0   0   0   0   0   D2 |
00166         //                  -                                             -     -                        -
00167         //
00168         //
00169         // We have something like
00170         //
00171         //  Bi^T D Bj
00172         //
00173         // These matrices have a particular pattern, especially the B-matrices:
00174         //
00175         //
00176         //                   |x 0 0|T
00177         //                   |0 x 0|
00178         //                   |0 0 x|
00179         //  | x 0 0 x 0 x|   |x x 0|
00180         //  | 0 x 0 x x 0| = |0 x x|              Bk(i,j)^T = Bk(j,i)
00181         //  | 0 0 x 0 x x|   |x 0 x|
00182         //
00183         // And the product has the pattern
00184         //
00185         //
00186         //                 | x x x 0 0 0 | |x 0 0|
00187         //                 | x x x 0 0 0 | |0 x 0|
00188         //                 | x x x 0 0 0 | |0 0 x|
00189         //  | x 0 0 x 0 x| | 0 0 0 x 0 0 | |x x 0|
00190         //  | 0 x 0 x x 0| | 0 0 0 0 x 0 | |0 x x|
00191         //  | 0 0 x 0 x x| | 0 0 0 0 0 x | |x 0 x|
00192         //
00193         //
00194         // If we compute T = Bi^T E, then we get a matrix with pattern
00195         //
00196         //  | x x x x 0 x|
00197         //  | x x x x x 0|
00198         //  | x x x 0 x x|
00199         //
00200         // Now computing  Kij = T Bj,
00201         //
00202         //                            |x 0 0 |
00203         //                            |0 x 0 |
00204         //                            |0 0 x |
00205         //  |x x x |   | x x x x 0 x| |x x 0 |
00206         //  |x x x | = | x x x x x 0| |0 x x |
00207         //  |x x x |   | x x x 0 x x| |x 0 x |
00208         //
00209         //
00210         // and exploiting that for isotropic materials we have,
00211         //
00212         //  D(0,0) = D(1,1) = D(2,2) = D0
00213         //  D(0,1) = D(1,2) = D(0,2) = D(1,0) = D(2,1) = D(2,0) = D1
00214         //  D(3,3) = D(4,4) = D(5,5) = D2
00215         //
00216         // Finally we do not need to store the B-matrices, observe that since
00217         //
00218         //  b(0) = invC(1,0);  b(1) = invC(1,1);  b(2) = invC(1,2);  b(3) = invC(1,3);
00219         //  c(0) = invC(2,0);  c(1) = invC(2,1);  c(2) = invC(2,2);  c(3) = invC(2,3);
00220         //  d(0) = invC(3,0);  d(1) = invC(3,1);  d(2) = invC(3,2);  d(3) = invC(3,3);
00221         //
00222         // and
00223         //
00224         //        | bi  0  0 |
00225         //        |  0 ci  0 |
00226         //   Bi = |  0  0 di |
00227         //        | ci bi  0 |
00228         //        |  0 di ci |
00229         //        | di  0 bi |
00230         //
00231         // therefore
00232         //
00233         //        | invC(1,i)    0         0      |
00234         //        |    0      invC(2,i)    0      |
00235         //   Bi = |    0         0      invC(3,i) |
00236         //        | invC(2,i) invC(1,i)    0      |
00237         //        |    0      invC(3,i) invC(2,i) |
00238         //        | invC(3,i)    0      invC(1,i) |
00239         //
00240         // Putting it all together and using matlab, we have,
00241         //
00242         //syms bi  ci  di real
00243         //Bi = [ bi  0  0 ;  0 ci  0 ;  0  0 di ; ci bi  0 ; 0 di ci ; di  0 bi ]
00244         //syms bj  cj  dj real
00245         //Bj = [ bj  0  0 ;  0 cj  0 ;  0  0 dj ; cj bj  0 ; 0 dj cj ; dj  0 bj ]
00246         //syms G H J real
00247         //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]
00248         //T = Bi'*E
00249         //Kij = T*Bj
00250         //syms invC1i invC2i invC3i real;
00251         //Kij = subs(Kij,bi,invC1i);
00252         //Kij = subs(Kij,ci,invC2i);
00253         //Kij = subs(Kij,di,invC3i);
00254         //syms invC1j invC2j invC3j real;
00255         //Kij = subs(Kij,bj,invC1j);
00256         //Kij = subs(Kij,cj,invC2j);
00257         //Kij = subs(Kij,dj,invC3j);
00258         //Kij = collect(Kij,G)
00259         //Kij = collect(Kij,H)
00260         //Kij = collect(Kij,J)
00261         //
00262         //  which yields
00263         //
00264         //             -                                                                                             -
00265         //             | D0 bi bj + D2 (ci cj + di dj)  D1 bi cj + D2 ci bj            D1 bi dj + D2 di bj           |
00266         //   Ke_ij = V | D1 ci bj + D2 bi cj            D0 ci cj + D2 (bi bj + di dj)  D1 ci dj + D2 di cj           |
00267         //             | D1 di bj + D2 bi dj            D1 di cj + D2 ci dj            D0 di dj + D2 (bi bj + ci cj) |
00268         //             -                                                                                             -
00269         real_type bi = Bi(0);
00270         real_type ci = Bi(1);
00271         real_type di = Bi(2);
00272 
00273         real_type bj = Bj(0);
00274         real_type cj = Bj(1);
00275         real_type dj = Bj(2);
00276 
00277         real_type D0 = D(0)*volume;
00278         real_type D1 = D(1)*volume;
00279         real_type D2 = D(2)*volume;
00280 
00281         Ke_ij(0,0) = D0 * bi * bj   +   D2 * (ci * cj + di * dj);
00282         Ke_ij(0,1) = D1 * bi * cj   +   D2 * (ci * bj);
00283         Ke_ij(0,2) = D1 * bi * dj   +   D2 * (di * bj);
00284 
00285         Ke_ij(1,0) = D1 * ci * bj   +   D2 * (bi * cj);
00286         Ke_ij(1,1) = D0 * ci * cj   +   D2 * (bi * bj + di * dj);
00287         Ke_ij(1,2) = D1 * ci * dj   +   D2 * (di * cj);
00288 
00289         Ke_ij(2,0) = D1 * di * bj   +   D2 * (bi * dj);
00290         Ke_ij(2,1) = D1 * di * cj   +   D2 * (ci * dj);
00291         Ke_ij(2,2) = D0 * di * dj   +   D2 * (ci * cj + bi * bj);
00292       }
00293 
00294     } // namespace detail
00295   } // namespace fem
00296 } // namespace OpenTissue
00297 
00298 //OPENTISSUE_DYNAMICS_FEM_FEM_COMPUTE_KE_H
00299 #endif
/home/hauberg/Dokumenter/Capture/humim-tracker-0.1/src/OpenTissue/OpenTissue/dynamics/fem/fem_compute_ke.h
