00001 #ifndef VIENNAFEM_TRANSFORM_TRIANGLE_HPP 00002 #define VIENNAFEM_TRANSFORM_TRIANGLE_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/triangle.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_triangle> 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 const & p0 = viennagrid::ncells<0>(cell)[0].point(); 00043 PointType const & p1 = viennagrid::ncells<0>(cell)[1].point(); 00044 PointType const & p2 = viennagrid::ncells<0>(cell)[2].point(); 00045 00046 //Step 1: store determinant: 00047 double det_dF_dt = 2.0 * viennagrid::spanned_volume(p0, p1, p2); 00048 00049 assert(det_dF_dt > 0); 00050 //std::cout << "dt_dx triangle!" << std::endl; 00051 00052 viennadata::access<det_dF_dt_key, numeric_type>()(cell) = det_dF_dt; 00053 00054 //Step 2: store partial derivatives: 00055 viennadata::access<dt_dx_key<0, 0>, numeric_type>()(cell) = ( p2[1] - p0[1]) / det_dF_dt; 00056 viennadata::access<dt_dx_key<0, 1>, numeric_type>()(cell) = - (p2[0] - p0[0]) / det_dF_dt; 00057 viennadata::access<dt_dx_key<1, 0>, numeric_type>()(cell) = - (p1[1] - p0[1]) / det_dF_dt; 00058 viennadata::access<dt_dx_key<1, 1>, numeric_type>()(cell) = ( p1[0] - p0[0]) / det_dF_dt; 00059 00060 } 00061 00062 }; 00063 00064 /* 00065 template <typename T_Configuration> 00066 struct dt_dx_handler<T_Configuration, TriangleTag, DtDxStoreDetOnly> 00067 { 00068 typedef typename T_Configuration::CoordType ScalarType; 00069 typedef typename DomainTypes<T_Configuration>::PointType PointType; 00070 typedef typename DomainTypes<T_Configuration>::VertexType VertexType; 00071 00072 public: 00073 dt_dx_handler() : Base() {}; 00074 dt_dx_handler( const dt_dx_handler & ddh ) : Base(ddh) {}; 00075 00076 //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. 00077 ScalarType get_dt_dx(int i, int j) const 00078 { 00079 PointType & p0 = Base::vertices_[0]->getPoint(); 00080 double ret = 0.0; 00081 00082 if (i == 0) 00083 { 00084 PointType & p2 = Base::vertices_[2]->getPoint(); 00085 if (j == 0) 00086 ret = (p2.get_y() - p0.get_y()) / det_dF_dt; 00087 else if (j == 1) 00088 ret = - (p2.get_x() - p0.get_x()) / det_dF_dt; 00089 else 00090 std::cerr << "ERROR: Accessing invalid elements of functional determinant!!"; 00091 } 00092 else if (i == 1) 00093 { 00094 PointType & p1 = Base::vertices_[1]->getPoint(); 00095 if (j == 0) 00096 ret = - (p1.get_y() - p0.get_y()) / det_dF_dt; 00097 else if (j == 1) 00098 ret = (p1.get_x() - p0.get_x()) / det_dF_dt; 00099 else 00100 std::cerr << "ERROR: Accessing invalid elements of functional determinant!!"; 00101 } 00102 00103 return ret; //dt_dx[2*i + j]; 00104 } 00105 00106 ScalarType get_det_dF_dt() const 00107 { 00108 return det_dF_dt; 00109 } 00110 00111 void update_dt_dx() 00112 { 00113 checkOrientation(); 00114 }; 00115 00116 void init_dt_dx() {} 00117 00118 protected: 00119 00120 void checkOrientation() 00121 { 00122 PointType & p0 = Base::vertices_[0]->getPoint(); 00123 PointType & p1 = Base::vertices_[1]->getPoint(); 00124 PointType & p2 = Base::vertices_[2]->getPoint(); 00125 00126 //volume: 00127 det_dF_dt = spannedVolume(p1 - p0, p2 - p0); 00128 00129 if (det_dF_dt < 0) 00130 { 00131 det_dF_dt = - det_dF_dt; 00132 VertexType *temp = Base::vertices_[0]; 00133 Base::vertices_[0] = Base::vertices_[1]; 00134 Base::vertices_[1] = temp; 00135 } 00136 else if (det_dF_dt == 0.0) 00137 std::cout << "ERROR: detected degenerated element!" << std::endl; 00138 } 00139 00140 void print(long indent) const 00141 { 00142 for (long i = 0; i<indent; ++i) 00143 std::cout << " "; 00144 std::cout << "* dt-dx-Handler: StoreDetOnly" << std::endl; 00145 Base::print(indent); 00146 } 00147 00148 private: 00149 ScalarType det_dF_dt; //determinant of Jacobian matrix 00150 }; 00151 00152 //save as much memory as possible: compute all values on access 00153 //however, in general the DtDxStoreDetOnly is much faster while having only moderate additional memory requirements. 00154 template <typename T_Configuration> 00155 struct dt_dx_handler<T_Configuration, TriangleTag, DtDxOnAccess> 00156 : public lower_level_holder<T_Configuration, TriangleTag, TriangleTag::TopoLevel - 1> 00157 { 00158 typedef lower_level_holder<T_Configuration, TriangleTag, TriangleTag::TopoLevel - 1> Base; 00159 typedef typename T_Configuration::CoordType ScalarType; 00160 typedef typename DomainTypes<T_Configuration>::PointType PointType; 00161 typedef typename DomainTypes<T_Configuration>::VertexType VertexType; 00162 00163 public: 00164 dt_dx_handler() : Base() {}; 00165 dt_dx_handler( const dt_dx_handler & ddh ) : Base(ddh) {}; 00166 00167 //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. 00168 ScalarType get_dt_dx(int i, int j) const 00169 { 00170 PointType & p0 = Base::vertices_[0]->getPoint(); 00171 PointType & p1 = Base::vertices_[1]->getPoint(); 00172 PointType & p2 = Base::vertices_[2]->getPoint(); 00173 double ret = 0.0; 00174 double det_dF_dt = spannedVolume(p1 - p0, p2 - p0); 00175 00176 if (i == 0) 00177 { 00178 if (j == 0) 00179 ret = (p2.get_y() - p0.get_y()) / det_dF_dt; 00180 else if (j == 1) 00181 ret = - (p2.get_x() - p0.get_x()) / det_dF_dt; 00182 else 00183 std::cerr << "ERROR: Accessing invalid elements of functional determinant!!"; 00184 } 00185 else if (i == 1) 00186 { 00187 if (j == 0) 00188 ret = - (p1.get_y() - p0.get_y()) / det_dF_dt; 00189 else if (j == 1) 00190 ret = (p1.get_x() - p0.get_x()) / det_dF_dt; 00191 else 00192 std::cerr << "ERROR: Accessing invalid elements of functional determinant!!"; 00193 } 00194 00195 return ret; //dt_dx[2*i + j]; 00196 } 00197 00198 ScalarType get_det_dF_dt() const 00199 { 00200 PointType & p0 = Base::vertices_[0]->getPoint(); 00201 PointType & p1 = Base::vertices_[1]->getPoint(); 00202 PointType & p2 = Base::vertices_[2]->getPoint(); 00203 return spannedVolume(p1 - p0, p2 - p0); 00204 } 00205 00206 void update_dt_dx() 00207 { 00208 checkOrientation(); 00209 }; 00210 00211 void init_dt_dx() {} 00212 00213 protected: 00214 00215 void checkOrientation() 00216 { 00217 PointType & p0 = Base::vertices_[0]->getPoint(); 00218 PointType & p1 = Base::vertices_[1]->getPoint(); 00219 PointType & p2 = Base::vertices_[2]->getPoint(); 00220 00221 //volume: 00222 double det_dF_dt = spannedVolume(p1 - p0, p2 - p0); 00223 00224 if (det_dF_dt < 0) 00225 { 00226 det_dF_dt = - det_dF_dt; 00227 VertexType *temp = Base::vertices_[0]; 00228 Base::vertices_[0] = Base::vertices_[1]; 00229 Base::vertices_[1] = temp; 00230 } 00231 } 00232 00233 void print(long indent = 0) const 00234 { 00235 for (long i = 0; i<indent; ++i) 00236 std::cout << " "; 00237 std::cout << "* dt-dx-Handler: OnAccess" << std::endl; 00238 Base::print(indent); 00239 } 00240 00241 }; 00242 00243 //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!! 00244 template <typename T_Configuration> 00245 struct dt_dx_handler<T_Configuration, TriangleTag, DtDxStoreStatically> 00246 : public lower_level_holder<T_Configuration, TriangleTag, TriangleTag::TopoLevel - 1> 00247 { 00248 typedef lower_level_holder<T_Configuration, TriangleTag, TriangleTag::TopoLevel - 1> Base; 00249 typedef typename T_Configuration::CoordType ScalarType; 00250 typedef typename DomainTypes<T_Configuration>::PointType PointType; 00251 typedef typename DomainTypes<T_Configuration>::VertexType VertexType; 00252 00253 public: 00254 dt_dx_handler() : Base() {}; 00255 dt_dx_handler( const dt_dx_handler & ddh ) : Base(ddh) {}; 00256 00257 //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. 00258 ScalarType get_dt_dx(int i, int j) const 00259 { 00260 return dt_dx[2*i + j]; 00261 } 00262 00263 ScalarType get_det_dF_dt() const 00264 { 00265 return det_dF_dt; 00266 } 00267 00268 void update_dt_dx() 00269 { 00270 checkOrientation(); 00271 }; 00272 00273 void init_dt_dx() 00274 { 00275 computeCellJacobian(); 00276 } 00277 00278 protected: 00279 00280 void computeCellJacobian() 00281 { 00282 PointType & p0 = Base::vertices_[0]->getPoint(); 00283 PointType & p1 = Base::vertices_[1]->getPoint(); 00284 PointType & p2 = Base::vertices_[2]->getPoint(); 00285 00286 det_dF_dt = spannedVolume(p1 - p0, p2 - p0); 00287 00288 if (det_dF_dt != 0) 00289 { 00290 //dt_1/dx_1 00291 dt_dx[0] = ( p2.get_y() - p0.get_y()) / det_dF_dt; 00292 //dt_1/dx_2 00293 dt_dx[1] = - (p2.get_x() - p0.get_x()) / det_dF_dt; 00294 00295 //dt_2/dx_1 00296 dt_dx[2] = - (p1.get_y() - p0.get_y()) / det_dF_dt; 00297 //dt_2/dx_2 00298 dt_dx[3] = ( p1.get_x() - p0.get_x()) / det_dF_dt; 00299 } 00300 00301 } //computeCellJacobian() 00302 00303 void checkOrientation() 00304 { 00305 PointType & p0 = Base::vertices_[0]->getPoint(); 00306 PointType & p1 = Base::vertices_[1]->getPoint(); 00307 PointType & p2 = Base::vertices_[2]->getPoint(); 00308 00309 //volume: 00310 ScalarType detdFdt = spannedVolume(p1 - p0, p2 - p0); 00311 00312 if (detdFdt < 0) 00313 { 00314 VertexType *temp = Base::vertices_[0]; 00315 Base::vertices_[0] = Base::vertices_[1]; 00316 Base::vertices_[1] = temp; 00317 } 00318 else if (detdFdt == 0.0) 00319 { 00320 std::cout << "ERROR: detected degenerated element!" << std::endl; 00321 Base::vertices_[0]->print(); 00322 Base::vertices_[1]->print(); 00323 Base::vertices_[2]->print(); 00324 std::cout << "Enter char to continue" << std::endl; 00325 char dummy; 00326 std::cin >> dummy; 00327 } 00328 } 00329 00330 void print(long indent) const 00331 { 00332 for (long i = 0; i<indent; ++i) 00333 std::cout << " "; 00334 std::cout << "* dt-dx-Handler: StoreStatically" << std::endl; 00335 Base::print(indent); 00336 } 00337 00338 private: 00339 static ScalarType dt_dx[4]; //inverse of Jacobian matrix from mapping 00340 static ScalarType det_dF_dt; //determinant of Jacobian matrix 00341 }; 00342 00343 //initialize static variables: 00344 template <typename T_Configuration> 00345 typename T_Configuration::CoordType 00346 dt_dx_handler<T_Configuration, TriangleTag, DtDxStoreStatically>::det_dF_dt = 0; 00347 00348 template <typename T_Configuration> 00349 typename T_Configuration::CoordType 00350 dt_dx_handler<T_Configuration, TriangleTag, DtDxStoreStatically>::dt_dx[] = {0, 0, 0, 0}; 00351 */ 00352 00353 } //namespace 00354 00355 #endif
1.7.1