00001 #ifndef VIENNAFEM_TRANSFORM_TETRAHEDRON_HPP 00002 #define VIENNAFEM_TRANSFORM_TETRAHEDRON_HPP 00003 00004 /* ========================================================================= 00005 Copyright (c) 2012, Institute for Microelectronics, 00006 Institute for Analysis and Scientific Computing, 00007 TU Wien. 00008 ----------------- 00009 ViennaFEM - The Vienna Finite Element Method Library 00010 ----------------- 00011 00012 Author: Karl Rupp rupp@iue.tuwien.ac.at 00013 00014 License: MIT (X11), see file LICENSE in the ViennaFEM base directory 00015 ============================================================================ */ 00016 00017 #include <iostream> 00018 #include "viennagrid/topology/tetrahedron.hpp" 00019 #include "viennagrid/algorithm/spanned_volume.hpp" 00020 #include "viennagrid/domain.hpp" 00021 #include "viennafem/forwards.h" 00022 00027 namespace viennafem 00028 { 00029 00030 //memory-intensive: Compute them once and store the computed values until next update 00031 template <> 00032 struct dt_dx_handler<viennafem::unit_tetrahedron> 00033 { 00034 public: 00035 00036 template <typename CellType> 00037 static void apply(CellType const & cell) 00038 { 00039 typedef typename CellType::config_type Config; 00040 typedef typename viennagrid::result_of::point<Config>::type PointType; 00041 00042 PointType p0 = viennagrid::ncells<0>(cell)[0].point(); 00043 PointType p1 = viennagrid::ncells<0>(cell)[1].point() - p0; 00044 PointType p2 = viennagrid::ncells<0>(cell)[2].point() - p0; 00045 PointType p3 = viennagrid::ncells<0>(cell)[3].point() - p0; 00046 00047 //Step 1: store determinant: 00048 numeric_type det_dF_dt = 6.0 * viennagrid::spanned_volume(viennagrid::ncells<0>(cell)[0].point(), 00049 viennagrid::ncells<0>(cell)[1].point(), 00050 viennagrid::ncells<0>(cell)[2].point(), 00051 viennagrid::ncells<0>(cell)[3].point()); 00052 00053 viennadata::access<det_dF_dt_key, numeric_type>()(cell) = det_dF_dt; 00054 00055 //Step 2: store partial derivatives: 00056 viennadata::access<dt_dx_key<0, 0>, numeric_type>()(cell) = ( + p2[1] * p3[2] - p2[2] * p3[1] ) / det_dF_dt; 00057 viennadata::access<dt_dx_key<0, 1>, numeric_type>()(cell) = ( - p2[0] * p3[2] + p2[2] * p3[0] ) / det_dF_dt; 00058 viennadata::access<dt_dx_key<0, 2>, numeric_type>()(cell) = ( + p2[0] * p3[1] - p2[1] * p3[0] ) / det_dF_dt; 00059 00060 viennadata::access<dt_dx_key<1, 0>, numeric_type>()(cell) = ( - p1[1] * p3[2] + p1[2] * p3[1] ) / det_dF_dt; 00061 viennadata::access<dt_dx_key<1, 1>, numeric_type>()(cell) = ( + p1[0] * p3[2] - p1[2] * p3[0] ) / det_dF_dt; 00062 viennadata::access<dt_dx_key<1, 2>, numeric_type>()(cell) = ( - p1[0] * p3[1] + p1[1] * p3[0] ) / det_dF_dt; 00063 00064 viennadata::access<dt_dx_key<2, 0>, numeric_type>()(cell) = ( + p1[1] * p2[2] - p1[2] * p2[1] ) / det_dF_dt; 00065 viennadata::access<dt_dx_key<2, 1>, numeric_type>()(cell) = ( - p1[0] * p2[2] + p1[2] * p2[0] ) / det_dF_dt; 00066 viennadata::access<dt_dx_key<2, 2>, numeric_type>()(cell) = ( + p1[0] * p2[1] - p1[1] * p2[0] ) / det_dF_dt; 00067 } 00068 00069 }; 00070 00071 00072 00073 /* Old code to follow... 00074 00075 00076 //memory-intensive: Compute them once and store the computed values until next update 00077 template <typename T_Configuration> 00078 struct dt_dx_handler<T_Configuration, viennagrid::tetrahedron_tag, DtDxStoreAll> 00079 { 00080 typedef typename T_Configuration::CoordType ScalarType; 00081 typedef typename viennagrid::DomainTypes<T_Configuration>::point_type PointType; 00082 typedef typename viennagrid::DomainTypes<T_Configuration>::vertex_type VertexType; 00083 00084 public: 00085 00086 //returns the element dt_i/dx_j of the functional determinant induced by the mapping to the reference element. i and j start at 0. 00087 ScalarType get_dt_dx(int i, int j) const 00088 { 00089 return dt_dx[3*i + j]; 00090 } 00091 00092 ScalarType get_det_dF_dt() const 00093 { 00094 return det_dF_dt; 00095 } 00096 00097 template <typename CellType> 00098 void update_dt_dx(CellType const & cell) 00099 { 00100 //checkOrientation(); 00101 00102 PointType p0 = cell.getPoint(0); 00103 PointType p1 = cell.getPoint(1) - p0; 00104 PointType p2 = cell.getPoint(2) - p0; 00105 PointType p3 = cell.getPoint(3) - p0; 00106 00107 //dt_1/dx_1 00108 dt_dx[0] = ( + p2[1] * p3[2] - p2[2] * p3[1] ) / det_dF_dt; 00109 //dt_1/dx_2 00110 dt_dx[1] = ( - p2[0] * p3[2] + p2[2] * p3[0] ) / det_dF_dt; 00111 //dt_1/dx_3 00112 dt_dx[2] = ( + p2[0] * p3[1] - p2[1] * p3[0] ) / det_dF_dt; 00113 00114 //dt_2/dx_1 00115 dt_dx[3] = ( - p1[1] * p3[2] + p1[2] * p3[1] ) / det_dF_dt; 00116 //dt_2/dx_2 00117 dt_dx[4] = ( + p1[0] * p3[2] - p1[2] * p3[0] ) / det_dF_dt; 00118 //dt_2/dx_3 00119 dt_dx[5] = ( - p1[0] * p3[1] + p1[1] * p3[0] ) / det_dF_dt; 00120 00121 //dt_3/dx_1 00122 dt_dx[6] = ( + p1[1] * p2[2] - p1[2] * p2[1] ) / det_dF_dt; 00123 //dt_3/dx_2 00124 dt_dx[7] = ( - p1[0] * p2[2] + p1[2] * p2[0] ) / det_dF_dt; 00125 //dt_3/dx_3 00126 dt_dx[8] = ( + p1[0] * p2[1] - p1[1] * p2[0] ) / det_dF_dt; 00127 00128 } 00129 00130 void init_dt_dx() {} 00131 00132 protected: 00133 00134 / * 00135 void checkOrientation() 00136 { 00137 PointType & p0 = cell.getPoint(0); 00138 PointType & p1 = cell.getPoint(1) - p0; 00139 PointType & p2 = cell.getPoint(2) - p0; 00140 PointType & p3 = cell.getPoint(3) - p0; 00141 00142 //volume: 00143 det_dF_dt = spannedVolume(p1 - p0, p2 - p0, p3 - p0); 00144 00145 if (det_dF_dt < 0) 00146 { 00147 det_dF_dt = - det_dF_dt; 00148 VertexType *temp = Base::vertices_[0]; 00149 Base::vertices_[0] = Base::vertices_[1]; 00150 Base::vertices_[1] = temp; 00151 } 00152 else if (det_dF_dt == 0.0) 00153 std::cout << "ERROR: detected degenerated element!" << std::endl; 00154 } * / 00155 00156 void print(long indent) const 00157 { 00158 for (long i = 0; i<indent; ++i) 00159 std::cout << " "; 00160 std::cout << "* dt-dx-Handler: StoreAll" << std::endl; 00161 } 00162 00163 private: 00164 ScalarType dt_dx[9]; //inverse of Jacobian matrix from mapping 00165 ScalarType det_dF_dt; //determinant of Jacobian matrix 00166 }; 00167 00168 //memory-computation-tradeoff: Store value of Jacobian only 00169 template <typename T_Configuration> 00170 struct dt_dx_handler<T_Configuration, viennagrid::tetrahedron_tag, DtDxStoreDetOnly> 00171 { 00172 typedef typename T_Configuration::numeric_type ScalarType; 00173 typedef typename viennagrid::DomainTypes<T_Configuration>::point_type PointType; 00174 typedef typename viennagrid::DomainTypes<T_Configuration>::vertex_type VertexType; 00175 00176 public: 00177 00178 //returns the element dt_i/dx_j of the functional determinant induced by the mapping to the reference element. i and j start at 0. 00179 template <typename CellType> 00180 ScalarType get_dt_dx(CellType const & cell, int i, int j) const 00181 { 00182 PointType p0 = cell.getPoint(0); 00183 PointType p1 = cell.getPoint(1) - p0; 00184 PointType p2 = cell.getPoint(2) - p0; 00185 PointType p3 = cell.getPoint(3) - p0; 00186 00187 double ret = 0.0; 00188 00189 //dt_1/dx_1 00190 if (i==0) 00191 { 00192 if (j==0) 00193 ret = ( + p2[1] * p3[2] - p2[2] * p3[1] ) / det_dF_dt; 00194 else if (j==1) 00195 ret = ( - p2[0] * p3[2] + p2[2] * p3[0] ) / det_dF_dt; 00196 else if (j==2) 00197 ret = ( + p2[0] * p3[1] - p2[1] * p3[0] ) / det_dF_dt; 00198 else 00199 std::cerr << "ERROR: Accessing invalid elements of functional determinant!!"; 00200 } 00201 else if (i==1) 00202 { 00203 if (j==0) 00204 ret = ( - p1[1] * p3[2] + p1[2] * p3[1] ) / det_dF_dt; 00205 else if (j==1) 00206 ret = ( + p1[0] * p3[2] - p1[2] * p3[0] ) / det_dF_dt; 00207 else if (j==2) 00208 ret = ( - p1[0] * p3[1] + p1[1] * p3[0] ) / det_dF_dt; 00209 else 00210 std::cerr << "ERROR: Accessing invalid elements of functional determinant!!"; 00211 } 00212 else if (i==2) 00213 { 00214 if (j==0) 00215 ret = ( + p1[1] * p2[2] - p1[2] * p2[1] ) / det_dF_dt; 00216 else if (j==1) 00217 ret = ( - p1[0] * p2[2] + p1[2] * p2[0] ) / det_dF_dt; 00218 else if (j==2) 00219 ret = ( + p1[0] * p2[1] - p1[1] * p2[0] ) / det_dF_dt; 00220 else 00221 std::cerr << "ERROR: Accessing invalid elements of functional determinant!!"; 00222 } 00223 else 00224 std::cerr << "ERROR: Accessing invalid elements of functional determinant!!"; 00225 00226 return ret; 00227 } 00228 00229 ScalarType get_det_dF_dt() const 00230 { 00231 return det_dF_dt; 00232 } 00233 00234 void update_dt_dx() 00235 { 00236 //checkOrientation(); 00237 } 00238 00239 void init_dt_dx() {} 00240 00241 protected: 00242 00243 / * 00244 void checkOrientation() 00245 { 00246 PointType & p0 = Base::vertices_[0]->getPoint(); 00247 PointType & p1 = Base::vertices_[1]->getPoint(); 00248 PointType & p2 = Base::vertices_[2]->getPoint(); 00249 PointType & p3 = Base::vertices_[3]->getPoint(); 00250 00251 //volume: 00252 det_dF_dt = spannedVolume(p1 - p0, p2 - p0, p3 - p0); 00253 00254 if (det_dF_dt < 0) 00255 { 00256 det_dF_dt = - det_dF_dt; 00257 VertexType *temp = Base::vertices_[0]; 00258 Base::vertices_[0] = Base::vertices_[1]; 00259 Base::vertices_[1] = temp; 00260 } 00261 else if (det_dF_dt == 0.0) 00262 std::cout << "ERROR: detected degenerated element!" << std::endl; 00263 } * / 00264 00265 void print(long indent) const 00266 { 00267 for (long i = 0; i<indent; ++i) 00268 std::cout << " "; 00269 std::cout << "* dt-dx-Handler: StoreDetOnly"<< std::endl; 00270 } 00271 00272 private: 00273 ScalarType det_dF_dt; //determinant of Jacobian matrix 00274 }; 00275 00276 //save as much memory as possible: compute all values on access 00277 //however, in general the DtDxStoreDetOnly is much faster while having only moderate additional memory requirements. 00278 template <typename T_Configuration> 00279 struct dt_dx_handler<T_Configuration, viennagrid::tetrahedron_tag, DtDxOnAccess> 00280 { 00281 typedef typename T_Configuration::numeric_type ScalarType; 00282 typedef typename viennagrid::DomainTypes<T_Configuration>::point_type PointType; 00283 typedef typename viennagrid::DomainTypes<T_Configuration>::vertex_type VertexType; 00284 00285 public: 00286 00287 //returns the element dt_i/dx_j of the functional determinant induced by the mapping to the reference element. i and j start at 0. 00288 template <typename CellType> 00289 ScalarType get_dt_dx(CellType const & cell, int i, int j) const 00290 { 00291 PointType p0 = cell.getPoint(0); 00292 PointType p1 = cell.getPoint(1) - p0; 00293 PointType p2 = cell.getPoint(2) - p0; 00294 PointType p3 = cell.getPoint(3) - p0; 00295 00296 double ret = 0.0; 00297 double det_dF_dt = get_det_dF_dt(); 00298 00299 //dt_1/dx_1 00300 if (i==0) 00301 { 00302 if (j==0) 00303 ret = ( + p2[1] * p3[2] - p2[2] * p3[1] ) / det_dF_dt; 00304 else if (j==1) 00305 ret = ( - p2[0] * p3[2] + p2[2] * p3[0] ) / det_dF_dt; 00306 else if (j==2) 00307 ret = ( + p2[0] * p3[1] - p2[1] * p3[0] ) / det_dF_dt; 00308 else 00309 std::cerr << "ERROR: Accessing invalid elements of functional determinant!!"; 00310 } 00311 else if (i==1) 00312 { 00313 if (j==0) 00314 ret = ( - p1[1] * p3[2] + p1[2] * p3[1] ) / det_dF_dt; 00315 else if (j==1) 00316 ret = ( + p1[0] * p3[2] - p1[2] * p3[0] ) / det_dF_dt; 00317 else if (j==2) 00318 ret = ( - p1[0] * p3[1] + p1[1] * p3[0] ) / det_dF_dt; 00319 else 00320 std::cerr << "ERROR: Accessing invalid elements of functional determinant!!"; 00321 } 00322 else if (i==2) 00323 { 00324 if (j==0) 00325 ret = ( + p1[1] * p2[2] - p1[2] * p2[1] ) / det_dF_dt; 00326 else if (j==1) 00327 ret = ( - p1[0] * p2[2] + p1[2] * p2[0] ) / det_dF_dt; 00328 else if (j==2) 00329 ret = ( + p1[0] * p2[1] - p1[1] * p2[0] ) / det_dF_dt; 00330 else 00331 std::cerr << "ERROR: Accessing invalid elements of functional determinant!!"; 00332 } 00333 else 00334 std::cerr << "ERROR: Accessing invalid elements of functional determinant!!"; 00335 00336 return ret; 00337 } 00338 00339 template <typename CellType> 00340 ScalarType get_det_dF_dt(CellType const & cell) const 00341 { 00342 PointType p0 = cell.getPoint(0); 00343 PointType p1 = cell.getPoint(1); 00344 PointType p2 = cell.getPoint(2); 00345 PointType p3 = cell.getPoint(3); 00346 00347 //volume: 00348 return spannedVolume(p1 - p0, p2 - p0, p3 - p0); 00349 } 00350 00351 protected: 00352 00353 / * 00354 void checkOrientation() 00355 { 00356 PointType & p0 = Base::vertices_[0]->getPoint(); 00357 PointType & p1 = Base::vertices_[1]->getPoint(); 00358 PointType & p2 = Base::vertices_[2]->getPoint(); 00359 PointType & p3 = Base::vertices_[3]->getPoint(); 00360 00361 //volume: 00362 double det_dF_dt = spannedVolume(p1 - p0, p2 - p0, p3 - p0); 00363 00364 if (det_dF_dt < 0) 00365 { 00366 det_dF_dt = - det_dF_dt; 00367 VertexType *temp = Base::vertices_[0]; 00368 Base::vertices_[0] = Base::vertices_[1]; 00369 Base::vertices_[1] = temp; 00370 } 00371 else if (det_dF_dt == 0.0) 00372 std::cout << "ERROR: detected degenerated element!" << std::endl; 00373 } * / 00374 00375 void print(long indent) const 00376 { 00377 for (long i = 0; i<indent; ++i) 00378 std::cout << " "; 00379 std::cout << "* dt-dx-Handler: OnAccess" << std::endl; 00380 } 00381 00382 }; 00383 00384 / * //fast and memory-saving: Compute an element's Jacobian as soon as the first entry is accessed. Use static memory, so that only one Jacobian set is stored at the same time. Must not be used in multi-threaded applications!! 00385 template <typename T_Configuration> 00386 struct dt_dx_handler<T_Configuration, TetrahedronTag, DtDxStoreStatically> 00387 { 00388 typedef typename T_Configuration::CoordType ScalarType; 00389 typedef typename DomainTypes<T_Configuration>::PointType PointType; 00390 typedef typename DomainTypes<T_Configuration>::VertexType VertexType; 00391 00392 public: 00393 dt_dx_handler() : Base() {}; 00394 dt_dx_handler( const dt_dx_handler & ddh ) : Base(ddh) {}; 00395 00396 //returns the element dt_i/dx_j of the functional determinant induced by the mapping to the reference element. i and j start at 0. 00397 ScalarType get_dt_dx(int i, int j) const 00398 { 00399 return dt_dx[3*i + j]; 00400 } 00401 00402 ScalarType get_det_dF_dt() const 00403 { 00404 return det_dF_dt; 00405 } 00406 00407 void update_dt_dx() 00408 { 00409 checkOrientation(); 00410 } 00411 00412 void print(long indent = 0) const 00413 { 00414 for (long i = 0; i<indent; ++i) 00415 std::cout << " "; 00416 std::cout << "* dt-dx-Handler: StoreStatically" << std::endl; 00417 Base::print(indent); 00418 } 00419 00420 void init_dt_dx() { computeCellJacobian(); } 00421 00422 protected: 00423 00424 void computeCellJacobian() const 00425 { 00426 PointType & p0 = Base::vertices_[0]->getPoint(); 00427 PointType p1 = Base::vertices_[1]->getPoint() - p0; 00428 PointType p2 = Base::vertices_[2]->getPoint() - p0; 00429 PointType p3 = Base::vertices_[3]->getPoint() - p0; 00430 00431 det_dF_dt = spannedVolume(p1, p2, p3); 00432 00433 //dt_1/dx_1 00434 dt_dx[0] = ( + p2[1] * p3[2] - p2[2] * p3[1] ) / det_dF_dt; 00435 //dt_1/dx_2 00436 dt_dx[1] = ( - p2[0] * p3[2] + p2[2] * p3[0] ) / det_dF_dt; 00437 //dt_1/dx_3 00438 dt_dx[2] = ( + p2[0] * p3[1] - p2[1] * p3[0] ) / det_dF_dt; 00439 00440 //dt_2/dx_1 00441 dt_dx[3] = ( - p1[1] * p3[2] + p1[2] * p3[1] ) / det_dF_dt; 00442 //dt_2/dx_2 00443 dt_dx[4] = ( + p1[0] * p3[2] - p1[2] * p3[0] ) / det_dF_dt; 00444 //dt_2/dx_3 00445 dt_dx[5] = ( - p1[0] * p3[1] + p1[1] * p3[0] ) / det_dF_dt; 00446 00447 //dt_3/dx_1 00448 dt_dx[6] = ( + p1[1] * p2[2] - p1[2] * p2[1] ) / det_dF_dt; 00449 //dt_3/dx_2 00450 dt_dx[7] = ( - p1[0] * p2[2] + p1[2] * p2[0] ) / det_dF_dt; 00451 //dt_3/dx_3 00452 dt_dx[8] = ( + p1[0] * p2[1] - p1[1] * p2[0] ) / det_dF_dt; 00453 00454 } //computeCellJacobian() 00455 00456 void checkOrientation() 00457 { 00458 PointType & p0 = Base::vertices_[0]->getPoint(); 00459 PointType & p1 = Base::vertices_[1]->getPoint(); 00460 PointType & p2 = Base::vertices_[2]->getPoint(); 00461 PointType & p3 = Base::vertices_[3]->getPoint(); 00462 00463 //volume: 00464 ScalarType detdFdt = spannedVolume(p1 - p0, p2 - p0, p3 - p0); 00465 00466 if (detdFdt < 0) 00467 { 00468 VertexType *temp = Base::vertices_[0]; 00469 Base::vertices_[0] = Base::vertices_[1]; 00470 Base::vertices_[1] = temp; 00471 } 00472 else if (detdFdt == 0.0) 00473 std::cout << "ERROR: detected degenerated element!" << std::endl; 00474 } 00475 00476 private: 00477 static ScalarType dt_dx[9]; //inverse of Jacobian matrix from mapping 00478 static ScalarType det_dF_dt; //determinant of Jacobian matrix 00479 }; 00480 00481 //initialize static variables: 00482 template <typename T_Configuration> 00483 typename T_Configuration::CoordType 00484 dt_dx_handler<T_Configuration, TetrahedronTag, DtDxStoreStatically>::det_dF_dt = 0; 00485 00486 template <typename T_Configuration> 00487 typename T_Configuration::CoordType 00488 dt_dx_handler<T_Configuration, TetrahedronTag, DtDxStoreStatically>::dt_dx[] = {0, 0, 0, 00489 0, 0, 0, 00490 0, 0, 0}; 00491 00492 00493 */ 00494 } 00495 #endif
1.7.1