00001 #ifndef OPENTISSUE_DYNAMICS_FEM_FEM_COMPUTE_B_H 00002 #define OPENTISSUE_DYNAMICS_FEM_FEM_COMPUTE_B_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 { 00018 00028 template<typename real_type, typename vector3_type> 00029 inline void compute_B( 00030 vector3_type const& e10, 00031 vector3_type const& e20, 00032 vector3_type const& e30, 00033 real_type volume, 00034 vector3_type * B 00035 ) 00036 { 00037 // 00038 // In summary we want to calculate 00039 // 00040 // B = S N 00041 // 00042 // where 00043 // 00044 // B = [ B0 B1 B2 B3 ], N = [ N0 N1 N2 N3 ] 00045 // 00046 // Here S is the operational matrix given by 00047 // - - 00048 // | d/dx 0 0 | 00049 // | 0 d/dy 0 | 00050 // | 0 0 d/dz | 00051 // S = | d/dy d/dx 0 | 00052 // | d/dz 0 d/dx | 00053 // | 0 d/dz d/dy | 00054 // - - 00055 // N_i's are called the shape functions and they are used to interpolate values. For 00056 // 3D linear elastostatics these can be written as 00057 // 00058 // - - 00059 // | w_i 0 0 | 00060 // N_i = | 0 w_i 0 | 00061 // | 0 0 w_i | 00062 // - - 00063 // 00064 // The w_i's are functions in the coordinates x,y, and z and will be derived 00065 // shortly below of here. 00066 // 00067 // The matrix B_i is given by: 00068 // 00069 // | b_i 0 0 | 00070 // | 0 c_i 0 | 00071 // B_i = | 0 0 d_i | 00072 // | c_i b_i 0 | 00073 // | d_i 0 b_i | 00074 // | 0 d_i c_i | 00075 // 00076 // The entries are stored as 00077 // 00078 // B[i] = [ b_i c_i d_i ] 00079 // 00080 // In the following we go into details of deriving formulas for the interpolating 00081 // functions N_i. From these formulas it follows how to compute the derivatives. 00082 // Finally we show a clever way of implementing the computaitons. 00083 // 00084 // Using a linear polynomial for interpolating a u(x,y,z)-value at a given 00085 // position x,y,z inside the tetrahedron 00086 // 00087 // u(x,y,z) = a_0 + a_1 x + a_2 y + a_3 z = P^T A 00088 // 00089 // where the a's are the polynomial coefficients (to be determed) and 00090 // 00091 // | 1 | | a0 | 00092 // P = | x | , and A = | a1 | 00093 // | y | | a2 | 00094 // | z | | a3 | 00095 // 00096 // Say we know the u-values, u0(x0,y0,z0),...,u3(x3,y2,z3), at the four tetrahedron corner 00097 // points, P0^T [x0 y0 z0],...,P3^T=[x3 y3 z3], then we can set up the linear system 00098 // 00099 // | u0 | | 1 x0 y0 z0 | | a0 | 00100 // | u1 | = | 1 x1 y1 z1 | | a1 | = C A 00101 // | u2 | | 1 x2 y2 z2 | | a2 | 00102 // | u3 | | 1 x3 y3 z3 | | a3 | 00103 // 00104 // The C matrix is invertible (as long as the four points are in general postion, meaning 00105 // that no point can be written as a linear combination of any of the three other points), 00106 // so we can solve the polynomial coefficients by 00107 // 00108 // | u0 | 00109 // A = C^-1 | u1 | 00110 // | u2 | 00111 // | u3 | 00112 // Substitution into the polynomial interpolation formula yields 00113 // 00114 // | u0 | 00115 // u(x,y,z) = P^T A = P^T C^{-1} | u1 | (*1) 00116 // | u2 | 00117 // | u3 | 00118 // 00119 // 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 00120 // 00121 // 00122 // Let us now try to look at the bary-centric coordinates, w0,...,w1 of (x,y,z), these 00123 // can also be used for interpolation 00124 // 00125 // | u0 | |u0| 00126 // u(z,y,z) = w0 u0 + w1 u1 + w2 u2 + w3 u3 = [w0 w1 w2 w3] | u1 | = W^T |u1| 00127 // | u2 | |u2| 00128 // | u3 | |u3| 00129 // 00130 // From the coordinates of the four points and the condition 1 = w0 + w1 + w2 + w3 we can set 00131 // up the linear system 00132 // 00133 // | 1 | | 1 1 1 1 | | w0 | 00134 // | x | = | x0 x1 x2 x3 | | w1 | 00135 // | y | | y0 y1 y2 y3 | | w2 | 00136 // | z | | z0 z1 z2 z3 | | w3 | 00137 // 00138 // 00139 // P = Q W 00140 // 00141 // Since the four points are in general postion, Q is invertible and we can solve for W 00142 // 00143 // W = Q^{-1} P 00144 // 00145 // Insertion into the barycentric interpolation formula yields 00146 // 00147 // |u0| |u0| |u0| 00148 // u(x,y,z) = W^T |u1| = ( Q^{-1} P )^T |u1| = P^T Q^{-T} |u1| (*2) 00149 // |u2| |u2| |u2| 00150 // |u3| |u3| |u3| 00151 // 00152 // Comparision with the polynomial interpolation derivation we see that C = Q^T, furthermore 00153 // we notice that w_i = N_i. So (not very surprinsingly) bary-centric interpolation is really 00154 // just linear polynomial interpolation. 00155 // 00156 // Notice that for the volume, V, of the tetrahedron we have the relation: 6 V = det( C ) = det ( Q^T ) 00157 // 00158 // When computing the stiffness matrix we are interested in derivatives of the N_i's with 00159 // respect to the x,y and z coordinates. 00160 // 00161 // 00162 // From (*1) and (*2) we see (recall we use zero-indexing) 00163 // 00164 // d N_i 00165 // ----- = C^{-1}_{1,i} = Q^{-T}_{1,i} (*3a) 00166 // d x 00167 // 00168 // d N_i 00169 // ----- = C^{-1}_{2,i} = Q^{-T}_{2,i} (*3b) 00170 // d y 00171 // 00172 // d N_i 00173 // ----- = C^{-1}_{3,i} = Q^{-T}_{3,i} (*3c) 00174 // d z 00175 // 00176 // Instead of actually computing the inverse of the 4x4 C or Q matrices a computational 00177 // more efficient solution can be derived, which only requires us to solve for a 3x3 system. 00178 // 00179 // By Cramers rule we have 00180 // 00181 // (-1)^{i+j} det( C_ji) 00182 // C^{-1}_{i,j} = ---------------------- 00183 // det(C) 00184 // 00185 // where det(C_ji) is the determinant of C with the j'th row and i'th column removed. 00186 // 00187 // Now defined 00188 // 00189 // e10 = P1 - P0 00190 // e20 = P2 - P0 00191 // e30 = P3 - P0 00192 // 00193 // and the matrix E 00194 // 00195 // | e10^T | | (x1-x0) (y1-y0) (z1-z0) | 00196 // E = | e20^T | = | (x2-x0) (y2-y0) (z2-z0) | 00197 // | e30^T | | (x3-x0) (y3-y0) (z3-z0) | 00198 // 00199 // Then 00200 // 00201 // det(C) = det(E) and C^{-1}_{i+1,j+1} = E^{-1}_{i,j} (*4) 00202 // 00203 // This can shown by straightforward computation, immediately is is verified that 00204 // 00205 // {-1}^((i+1)+(j+1)) = {-1}^(i+j) 00206 // 00207 // So we have to show 00208 // 00209 // det( C_(i+1,j+1)) det( E_(ij)) 00210 // --------------------- = -------------- 00211 // det(C) det(E) 00212 // 00213 // assume det(C)= det(E) (left for reader as exercise:-) then inorder to prove 00214 // second half of (*4) we have to show 00215 // 00216 // det( C_(i+1,j+1)) = det( E_(ij)) 00217 // 00218 // A total of 9 cases exist, for here we will show a single case and leave 00219 // the remaining cases for the reader, use i=1 and j= 0, implying 00220 // 00221 // det( C_(12)) = det( E_(01)) 00222 // 00223 // So 00224 // - - 00225 // | 1 x0 z0 | - - 00226 // det | 1 x2 z2 | = det | (x2-x0) (z2-z0)| 00227 // | 1 x3 z3 | | (x3-x0) (z3-z0)| 00228 // - - - - 00229 // 00230 // Which is trivially true. So we have 00231 // 00232 // | . ... | | . ... | 00233 // C^{-1} = | . E^{-1} | ; Q^{-T} = | . E^{-T} | 00234 // 00235 // As can be seen from (*3) we do not have all the values needed since 00236 // the first columns of C^{-1} and Q^{-T} are missing. To remedy this problem we can use 00237 // the condition of the bary-centric coordinates 00238 // 00239 // w0 = 1 - w1 - w2 - w3 00240 // 00241 // Taking the derivate wrt. x,y, and z yields 00242 // 00243 // d N_0 00244 // ----- = 1 - E^{-1}_01 - E^{-1}_02 - E^{-1}_03 (*5a) 00245 // d x 00246 // 00247 // d N_0 00248 // ----- = 1 - E^{-1}_11 - E^{-1}_12 - E^{-1}_13 (*5b) 00249 // d y 00250 // 00251 // d N_0 00252 // ----- = 1 - E^{-1}_21 - E^{-1}_22 - E^{-1}_23 (*5c) 00253 // d z 00254 // 00255 // 00256 // and for i>0 we have 00257 // 00258 // d N_i 00259 // ----- = E^{-1}_0i (*6a) 00260 // d x 00261 // 00262 // d N_i 00263 // ----- = E^{-1}_1i (*6b) 00264 // d y 00265 // 00266 // 00267 // d N_i 00268 // ----- = E^{-1}_2i (*6c) 00269 // d z 00270 // 00271 // The notation b_i = d N_i / dx, c_i = d N_i / dy, and d_i = d N_i / dz is used. Furthermore 00272 // all the derivatives are returned as four B-vectors, where 00273 // 00274 // B[i] = [ b_i c_i d_i] 00275 // 00276 // 00277 // 00278 // The above may seem as a mathematical trick, so let us try to derive our 00279 // equations a litlle differently, but before doing so we must first develop 00280 // some equations for the volume of a tetrahedron. 00281 // 00282 // | e10^T | 00283 // 6 V = det(E) = | e20^T | = e10 \cdot (e20 \times e30) 00284 // | e30^T | 00285 // 00286 // where e10 = p1-p0, e20 = p2-p0 and so on. 00287 // In general given the right-hand order i,j,k, and m of the nodes we 00288 // write vol(i,j,k,m) = eji \cdot (eki \times emi) / 6 00289 // 00290 // Barycentric coordinates are infact the volume weighted weights coresponding to 00291 // the four tetrahedra lying inside the tetrahedron and having apex at the 00292 // point p=[x,y,z]^T and bases equal to the 4 triangular faces of the enclosing 00293 // tetrahedron. To realize this we will examine bary-centric coordinates 00294 // in the 2D case, that is the case of the Triangle. In 2D area corresponds 00295 // to the volume, so given a trianle and a point p inside it 00296 // 00297 // + 2 00298 // /| 00299 // / | 00300 // / | 00301 // / | 00302 // / + p| 00303 // / | 00304 // / | 00305 // 0 +-------+ 1 00306 // 00307 // Let the area weighted weights be defined as 00308 // 00309 // w_0 = area(1,2,p) / area(0,1,2) 00310 // w_1 = area(0,P,2) / area(0,1,2) = area(2,0,P) / area(0,1,2) 00311 // w_2 = area(0,1,p) / area(0,1,2) 00312 // 00313 // It is intuitive to see that if p moves towards the i'th corner point then 00314 // the nominator of w_i goes towards area(0,1,2). This means that w_i -> 1 00315 // while w_j ->0 for j\neq i. 00316 // 00317 // Having introduced the barycentric coordinates as the area weighted weights 00318 // it is quite intuitive to see that we also must have 00319 // 00320 // 1 = \sum_i w_i 00321 // 00322 // In 3D the barycentric coordinates are the volume weighted weights. That is 00323 // given the right handed order 0,1,2 and 3 of the corner points, we have by 00324 // analogy to the 2D case 00325 // 00326 // vol(0,1,2,p) 00327 // w_3 = ------------ = (p-p0) \codt ( e10 \times \e20 ) / 6V 00328 // V 00329 // 00330 // vol(0,1,3,p) 00331 // w_2 = ------------ = (p-p0) \codt ( e10 \times \e30 ) / 6V 00332 // V 00333 // 00334 // vol(0,2,3,p) 00335 // w_1 = ------------ = (p-p0) \codt ( e20 \times \e30 ) / 6V 00336 // V 00337 // 00338 // vol(1,3,2,p) 00339 // w_0 = ------------ = (p-p0) \codt ( e31 \times \e21 ) / 6V 00340 // V 00341 // 00342 // But since we allready know w_3, w_2 and w_1 it is faster to compute w_0 as 00343 // 00344 // w_0 = 1 - w_1 - w_2 - w_3 00345 // 00346 // Recal from previous that we seek the derivatives of w_i's wrt. the coordinaes when 00347 // we compute the B functions. Let us write the [x y z] coordinates as [x_0 x_1 x_2] then 00348 // 00349 // d w_3 00350 // ------- = ( e10 \times \e20 )_{x_j} / 6V 00351 // d x_j 00352 // 00353 // That is the j'th coordinate of the cross product divided by 6 times the volume. Similar for 00354 // 00355 // d w_2 00356 // ------- = ( e10 \times \e30 )_{x_j} / 6V 00357 // d x_j 00358 // 00359 // d w_1 00360 // ------- = ( e20 \times \e30 )_{x_j} / 6V 00361 // d x_j 00362 // 00363 // and finally 00364 // 00365 // d w_0 d w_i 00366 // ------- = 1 - sum_i>0 ------- 00367 // d x_j d x_j 00368 // 00369 // Our formulas may seem differnt from those in (*6), however 00370 // 00371 // d w_k 00372 // ------- = ( eji \times \emi )_{x_j} / 6V = E_{-1}_{k,j} 00373 // d x_j 00374 // 00375 // That is the adjoint of E_(j,k) is given by ( eji \times \emi )_{x_j} / 6. 00376 // 00377 // 00378 // 00379 00380 real_type div6V = 1/(6.0*volume); 00381 00382 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) 00383 B[2](0) = (e10(1) * e30(2) - e10(2) * e30(1)) * div6V; // b1 = det(E_12) 00384 B[3](0) = (e10(2) * e20(1) - e10(1) * e20(2)) * div6V; // b2 = -det(E_13) 00385 B[0](0) = -B[1](0) - B[2](0) - B[3](0); // b3 = -b0 - b1 - b2 00386 00387 B[1](1) = (e20(0) * e30(2) - e20(2) * e30(0)) * div6V; // c0 = det(E_21) 00388 B[2](1) = (e10(2) * e30(0) - e10(0) * e30(2)) * div6V; // c1 = -det(E_22) 00389 B[3](1) = (e10(0) * e20(2) - e10(2) * e20(0)) * div6V; // c2 = det(E_23) 00390 B[0](1) = -B[1](1) - B[2](1) - B[3](1); // c3 = -c0 - c1 - c2 00391 00392 B[1](2) = (e20(1) * e30(0) - e20(0) * e30(1)) * div6V; // d0 = -det(E_31) 00393 B[2](2) = (e10(0) * e30(1) - e10(1) * e30(0)) * div6V; // d1 = det(E_32) 00394 B[3](2) = (e10(1) * e20(0) - e10(0) * e20(1)) * div6V; // d2 = -det(E_33) 00395 B[0](2) = -B[1](2) - B[2](2) - B[3](2); // d3 = -d0 - d1 - d2 00396 } 00397 00398 } // namespace detail 00399 } // namespace fem 00400 } // namespace OpenTissue 00401 00402 //OPENTISSUE_DYNAMICS_FEM_FEM_COMPUTE_B_H 00403 #endif
1.7.1