h1_fe_transformation.C
Go to the documentation of this file.00001 // The libMesh Finite Element Library. 00002 // Copyright (C) 2002-2012 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner 00003 00004 // This library is free software; you can redistribute it and/or 00005 // modify it under the terms of the GNU Lesser General Public 00006 // License as published by the Free Software Foundation; either 00007 // version 2.1 of the License, or (at your option) any later version. 00008 00009 // This library is distributed in the hope that it will be useful, 00010 // but WITHOUT ANY WARRANTY; without even the implied warranty of 00011 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 00012 // Lesser General Public License for more details. 00013 00014 // You should have received a copy of the GNU Lesser General Public 00015 // License along with this library; if not, write to the Free Software 00016 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA 00017 00018 #include "libmesh/fe_interface.h" 00019 #include "libmesh/h1_fe_transformation.h" 00020 #include "libmesh/tensor_value.h" 00021 00022 namespace libMesh 00023 { 00024 template< typename OutputShape > 00025 void H1FETransformation<OutputShape>::map_phi( const unsigned int dim, 00026 const Elem* const elem, 00027 const std::vector<Point>& qp, 00028 const FEGenericBase<OutputShape>& fe, 00029 std::vector<std::vector<OutputShape> >& phi ) const 00030 { 00031 switch(dim) 00032 { 00033 case 0: 00034 { 00035 for (unsigned int i=0; i<phi.size(); i++) 00036 { 00037 libmesh_assert_equal_to ( qp.size(), phi[i].size() ); 00038 for (unsigned int p=0; p<phi[i].size(); p++) 00039 FEInterface::shape<OutputShape>(0, fe.get_fe_type(), elem, i, qp[p], phi[i][p]); 00040 } 00041 break; 00042 } 00043 case 1: 00044 { 00045 for (unsigned int i=0; i<phi.size(); i++) 00046 { 00047 libmesh_assert_equal_to ( qp.size(), phi[i].size() ); 00048 for (unsigned int p=0; p<phi[i].size(); p++) 00049 FEInterface::shape<OutputShape>(1, fe.get_fe_type(), elem, i, qp[p], phi[i][p]); 00050 } 00051 break; 00052 } 00053 case 2: 00054 { 00055 for (unsigned int i=0; i<phi.size(); i++) 00056 { 00057 libmesh_assert_equal_to ( qp.size(), phi[i].size() ); 00058 for (unsigned int p=0; p<phi[i].size(); p++) 00059 FEInterface::shape<OutputShape>(2, fe.get_fe_type(), elem, i, qp[p], phi[i][p]); 00060 } 00061 break; 00062 } 00063 case 3: 00064 { 00065 for (unsigned int i=0; i<phi.size(); i++) 00066 { 00067 libmesh_assert_equal_to ( qp.size(), phi[i].size() ); 00068 for (unsigned int p=0; p<phi[i].size(); p++) 00069 FEInterface::shape<OutputShape>(3, fe.get_fe_type(), elem, i, qp[p], phi[i][p]); 00070 } 00071 break; 00072 } 00073 default: 00074 libmesh_error(); 00075 } 00076 00077 return; 00078 } 00079 00080 00081 template< typename OutputShape > 00082 void H1FETransformation<OutputShape>::map_dphi( const unsigned int dim, 00083 const Elem* const, 00084 const std::vector<Point>&, 00085 const FEGenericBase<OutputShape>& fe, 00086 std::vector<std::vector<typename FEGenericBase<OutputShape>::OutputGradient> >& dphi, 00087 std::vector<std::vector<OutputShape> >& dphidx, 00088 std::vector<std::vector<OutputShape> >& dphidy, 00089 std::vector<std::vector<OutputShape> >& dphidz ) const 00090 { 00091 switch(dim) 00092 { 00093 case 0: // No derivatives in 0D 00094 { 00095 for (unsigned int i=0; i<dphi.size(); i++) 00096 for (unsigned int p=0; p<dphi[i].size(); p++) 00097 { 00098 dphi[i][p] = 0.; 00099 } 00100 break; 00101 } 00102 00103 case 1: 00104 { 00105 const std::vector<std::vector<OutputShape> >& dphidxi = fe.get_dphidxi(); 00106 00107 const std::vector<Real>& dxidx_map = fe.get_fe_map().get_dxidx(); 00108 #if LIBMESH_DIM>1 00109 const std::vector<Real>& dxidy_map = fe.get_fe_map().get_dxidy(); 00110 #endif 00111 #if LIBMESH_DIM>2 00112 const std::vector<Real>& dxidz_map = fe.get_fe_map().get_dxidz(); 00113 #endif 00114 00115 for (unsigned int i=0; i<dphi.size(); i++) 00116 for (unsigned int p=0; p<dphi[i].size(); p++) 00117 { 00118 // dphi/dx = (dphi/dxi)*(dxi/dx) 00119 dphi[i][p].slice(0) = dphidx[i][p] = dphidxi[i][p]*dxidx_map[p]; 00120 00121 #if LIBMESH_DIM>1 00122 dphi[i][p].slice(1) = dphidy[i][p] = dphidxi[i][p]*dxidy_map[p]; 00123 #endif 00124 #if LIBMESH_DIM>2 00125 dphi[i][p].slice(2) = dphidz[i][p] = dphidxi[i][p]*dxidz_map[p]; 00126 #endif 00127 } 00128 00129 break; 00130 } 00131 00132 case 2: 00133 { 00134 const std::vector<std::vector<OutputShape> >& dphidxi = fe.get_dphidxi(); 00135 const std::vector<std::vector<OutputShape> >& dphideta = fe.get_dphideta(); 00136 00137 const std::vector<Real>& dxidx_map = fe.get_fe_map().get_dxidx(); 00138 const std::vector<Real>& dxidy_map = fe.get_fe_map().get_dxidy(); 00139 #if LIBMESH_DIM > 2 00140 const std::vector<Real>& dxidz_map = fe.get_fe_map().get_dxidz(); 00141 #endif 00142 00143 const std::vector<Real>& detadx_map = fe.get_fe_map().get_detadx(); 00144 const std::vector<Real>& detady_map = fe.get_fe_map().get_detady(); 00145 #if LIBMESH_DIM > 2 00146 const std::vector<Real>& detadz_map = fe.get_fe_map().get_detadz(); 00147 #endif 00148 00149 for (unsigned int i=0; i<dphi.size(); i++) 00150 for (unsigned int p=0; p<dphi[i].size(); p++) 00151 { 00152 // dphi/dx = (dphi/dxi)*(dxi/dx) + (dphi/deta)*(deta/dx) 00153 dphi[i][p].slice(0) = dphidx[i][p] = (dphidxi[i][p]*dxidx_map[p] + 00154 dphideta[i][p]*detadx_map[p]); 00155 00156 // dphi/dy = (dphi/dxi)*(dxi/dy) + (dphi/deta)*(deta/dy) 00157 dphi[i][p].slice(1) = dphidy[i][p] = (dphidxi[i][p]*dxidy_map[p] + 00158 dphideta[i][p]*detady_map[p]); 00159 00160 #if LIBMESH_DIM > 2 00161 // dphi/dz = (dphi/dxi)*(dxi/dz) + (dphi/deta)*(deta/dz) 00162 dphi[i][p].slice(2) = dphidz[i][p] = (dphidxi[i][p]*dxidz_map[p] + 00163 dphideta[i][p]*detadz_map[p]); 00164 #endif 00165 } 00166 00167 break; 00168 } 00169 00170 case 3: 00171 { 00172 const std::vector<std::vector<OutputShape> >& dphidxi = fe.get_dphidxi(); 00173 const std::vector<std::vector<OutputShape> >& dphideta = fe.get_dphideta(); 00174 const std::vector<std::vector<OutputShape> >& dphidzeta = fe.get_dphidzeta(); 00175 00176 const std::vector<Real>& dxidx_map = fe.get_fe_map().get_dxidx(); 00177 const std::vector<Real>& dxidy_map = fe.get_fe_map().get_dxidy(); 00178 const std::vector<Real>& dxidz_map = fe.get_fe_map().get_dxidz(); 00179 00180 const std::vector<Real>& detadx_map = fe.get_fe_map().get_detadx(); 00181 const std::vector<Real>& detady_map = fe.get_fe_map().get_detady(); 00182 const std::vector<Real>& detadz_map = fe.get_fe_map().get_detadz(); 00183 00184 const std::vector<Real>& dzetadx_map = fe.get_fe_map().get_dzetadx(); 00185 const std::vector<Real>& dzetady_map = fe.get_fe_map().get_dzetady(); 00186 const std::vector<Real>& dzetadz_map = fe.get_fe_map().get_dzetadz(); 00187 00188 for (unsigned int i=0; i<dphi.size(); i++) 00189 for (unsigned int p=0; p<dphi[i].size(); p++) 00190 { 00191 // dphi/dx = (dphi/dxi)*(dxi/dx) + (dphi/deta)*(deta/dx) + (dphi/dzeta)*(dzeta/dx); 00192 dphi[i][p].slice(0) = dphidx[i][p] = (dphidxi[i][p]*dxidx_map[p] + 00193 dphideta[i][p]*detadx_map[p] + 00194 dphidzeta[i][p]*dzetadx_map[p]); 00195 00196 // dphi/dy = (dphi/dxi)*(dxi/dy) + (dphi/deta)*(deta/dy) + (dphi/dzeta)*(dzeta/dy); 00197 dphi[i][p].slice(1) = dphidy[i][p] = (dphidxi[i][p]*dxidy_map[p] + 00198 dphideta[i][p]*detady_map[p] + 00199 dphidzeta[i][p]*dzetady_map[p]); 00200 00201 // dphi/dz = (dphi/dxi)*(dxi/dz) + (dphi/deta)*(deta/dz) + (dphi/dzeta)*(dzeta/dz); 00202 dphi[i][p].slice(2) = dphidz[i][p] = (dphidxi[i][p]*dxidz_map[p] + 00203 dphideta[i][p]*detadz_map[p] + 00204 dphidzeta[i][p]*dzetadz_map[p]); 00205 } 00206 break; 00207 } 00208 00209 default: 00210 libmesh_error(); 00211 00212 } // switch(dim) 00213 00214 return; 00215 } 00216 00217 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES 00218 template< typename OutputShape > 00219 void H1FETransformation<OutputShape>::map_d2phi( const unsigned int dim, 00220 const Elem* const elem, 00221 const std::vector<Point>&, 00222 const FEGenericBase<OutputShape>& fe, 00223 std::vector<std::vector<typename FEGenericBase<OutputShape>::OutputTensor> >& d2phi, 00224 std::vector<std::vector<OutputShape> >& d2phidx2, 00225 std::vector<std::vector<OutputShape> >& d2phidxdy, 00226 std::vector<std::vector<OutputShape> >& d2phidxdz, 00227 std::vector<std::vector<OutputShape> >& d2phidy2, 00228 std::vector<std::vector<OutputShape> >& d2phidydz, 00229 std::vector<std::vector<OutputShape> >& d2phidz2 ) const 00230 { 00231 libmesh_do_once( 00232 if (!elem->has_affine_map()) 00233 { 00234 libMesh::err << "WARNING: Second derivatives are not currently " 00235 << "correctly calculated on non-affine elements!" 00236 << std::endl; 00237 } 00238 ); 00239 00240 00241 switch(dim) 00242 { 00243 case 0: // No derivatives in 0D 00244 { 00245 for (unsigned int i=0; i<d2phi.size(); i++) 00246 for (unsigned int p=0; p<d2phi[i].size(); p++) 00247 { 00248 d2phi[i][p] = 0.; 00249 } 00250 00251 break; 00252 } 00253 00254 case 1: 00255 { 00256 const std::vector<std::vector<OutputShape> >& d2phidxi2 = fe.get_d2phidxi2(); 00257 00258 const std::vector<Real>& dxidx_map = fe.get_fe_map().get_dxidx(); 00259 #if LIBMESH_DIM>1 00260 const std::vector<Real>& dxidy_map = fe.get_fe_map().get_dxidy(); 00261 #endif 00262 #if LIBMESH_DIM>2 00263 const std::vector<Real>& dxidz_map = fe.get_fe_map().get_dxidz(); 00264 #endif 00265 00266 for (unsigned int i=0; i<d2phi.size(); i++) 00267 for (unsigned int p=0; p<d2phi[i].size(); p++) 00268 { 00269 d2phi[i][p].slice(0).slice(0) = d2phidx2[i][p] = 00270 d2phidxi2[i][p]*dxidx_map[p]*dxidx_map[p]; 00271 #if LIBMESH_DIM>1 00272 d2phi[i][p].slice(0).slice(1) = 00273 d2phi[i][p].slice(1).slice(0) = d2phidxdy[i][p] = 00274 d2phidxi2[i][p]*dxidx_map[p]*dxidy_map[p]; 00275 00276 d2phi[i][p].slice(1).slice(1) = d2phidy2[i][p] = 00277 d2phidxi2[i][p]*dxidy_map[p]*dxidy_map[p]; 00278 #endif 00279 #if LIBMESH_DIM>2 00280 d2phi[i][p].slice(0).slice(2) = 00281 d2phi[i][p].slice(2).slice(0) = d2phidxdz[i][p] = 00282 d2phidxi2[i][p]*dxidx_map[p]*dxidz_map[p]; 00283 00284 d2phi[i][p].slice(1).slice(2) = 00285 d2phi[i][p].slice(2).slice(1) = d2phidydz[i][p] = 00286 d2phidxi2[i][p]*dxidy_map[p]*dxidz_map[p]; 00287 00288 d2phi[i][p].slice(2).slice(2) = d2phidz2[i][p] = 00289 d2phidxi2[i][p]*dxidz_map[p]*dxidz_map[p]; 00290 #endif 00291 } 00292 break; 00293 } 00294 00295 case 2: 00296 { 00297 const std::vector<std::vector<OutputShape> >& d2phidxi2 = fe.get_d2phidxi2(); 00298 const std::vector<std::vector<OutputShape> >& d2phidxideta = fe.get_d2phidxideta(); 00299 const std::vector<std::vector<OutputShape> >& d2phideta2 = fe.get_d2phideta2(); 00300 00301 const std::vector<Real>& dxidx_map = fe.get_fe_map().get_dxidx(); 00302 const std::vector<Real>& dxidy_map = fe.get_fe_map().get_dxidy(); 00303 #if LIBMESH_DIM > 2 00304 const std::vector<Real>& dxidz_map = fe.get_fe_map().get_dxidz(); 00305 #endif 00306 00307 const std::vector<Real>& detadx_map = fe.get_fe_map().get_detadx(); 00308 const std::vector<Real>& detady_map = fe.get_fe_map().get_detady(); 00309 #if LIBMESH_DIM > 2 00310 const std::vector<Real>& detadz_map = fe.get_fe_map().get_detadz(); 00311 #endif 00312 00313 for (unsigned int i=0; i<d2phi.size(); i++) 00314 for (unsigned int p=0; p<d2phi[i].size(); p++) 00315 { 00316 d2phi[i][p].slice(0).slice(0) = d2phidx2[i][p] = 00317 d2phidxi2[i][p]*dxidx_map[p]*dxidx_map[p] + 00318 2*d2phidxideta[i][p]*dxidx_map[p]*detadx_map[p] + 00319 d2phideta2[i][p]*detadx_map[p]*detadx_map[p]; 00320 00321 d2phi[i][p].slice(0).slice(1) = 00322 d2phi[i][p].slice(1).slice(0) = d2phidxdy[i][p] = 00323 d2phidxi2[i][p]*dxidx_map[p]*dxidy_map[p] + 00324 d2phidxideta[i][p]*dxidx_map[p]*detady_map[p] + 00325 d2phideta2[i][p]*detadx_map[p]*detady_map[p] + 00326 d2phidxideta[i][p]*detadx_map[p]*dxidy_map[p]; 00327 00328 d2phi[i][p].slice(1).slice(1) = d2phidy2[i][p] = 00329 d2phidxi2[i][p]*dxidy_map[p]*dxidy_map[p] + 00330 2*d2phidxideta[i][p]*dxidy_map[p]*detady_map[p] + 00331 d2phideta2[i][p]*detady_map[p]*detady_map[p]; 00332 00333 #if LIBMESH_DIM > 2 00334 d2phi[i][p].slice(0).slice(2) = 00335 d2phi[i][p].slice(2).slice(0) = d2phidxdz[i][p] = 00336 d2phidxi2[i][p]*dxidx_map[p]*dxidz_map[p] + 00337 d2phidxideta[i][p]*dxidx_map[p]*detadz_map[p] + 00338 d2phideta2[i][p]*detadx_map[p]*detadz_map[p] + 00339 d2phidxideta[i][p]*detadx_map[p]*dxidz_map[p]; 00340 00341 d2phi[i][p].slice(1).slice(2) = 00342 d2phi[i][p].slice(2).slice(1) = d2phidydz[i][p] = 00343 d2phidxi2[i][p]*dxidy_map[p]*dxidz_map[p] + 00344 d2phidxideta[i][p]*dxidy_map[p]*detadz_map[p] + 00345 d2phideta2[i][p]*detady_map[p]*detadz_map[p] + 00346 d2phidxideta[i][p]*detady_map[p]*dxidz_map[p]; 00347 00348 d2phi[i][p].slice(2).slice(2) = d2phidz2[i][p] = 00349 d2phidxi2[i][p]*dxidz_map[p]*dxidz_map[p] + 00350 2*d2phidxideta[i][p]*dxidz_map[p]*detadz_map[p] + 00351 d2phideta2[i][p]*detadz_map[p]*detadz_map[p]; 00352 #endif 00353 } 00354 00355 break; 00356 } 00357 00358 case 3: 00359 { 00360 const std::vector<std::vector<OutputShape> >& d2phidxi2 = fe.get_d2phidxi2(); 00361 const std::vector<std::vector<OutputShape> >& d2phidxideta = fe.get_d2phidxideta(); 00362 const std::vector<std::vector<OutputShape> >& d2phideta2 = fe.get_d2phideta2(); 00363 const std::vector<std::vector<OutputShape> >& d2phidxidzeta = fe.get_d2phidxidzeta(); 00364 const std::vector<std::vector<OutputShape> >& d2phidetadzeta = fe.get_d2phidetadzeta(); 00365 const std::vector<std::vector<OutputShape> >& d2phidzeta2 = fe.get_d2phidzeta2(); 00366 00367 const std::vector<Real>& dxidx_map = fe.get_fe_map().get_dxidx(); 00368 const std::vector<Real>& dxidy_map = fe.get_fe_map().get_dxidy(); 00369 const std::vector<Real>& dxidz_map = fe.get_fe_map().get_dxidz(); 00370 00371 const std::vector<Real>& detadx_map = fe.get_fe_map().get_detadx(); 00372 const std::vector<Real>& detady_map = fe.get_fe_map().get_detady(); 00373 const std::vector<Real>& detadz_map = fe.get_fe_map().get_detadz(); 00374 00375 const std::vector<Real>& dzetadx_map = fe.get_fe_map().get_dzetadx(); 00376 const std::vector<Real>& dzetady_map = fe.get_fe_map().get_dzetady(); 00377 const std::vector<Real>& dzetadz_map = fe.get_fe_map().get_dzetadz(); 00378 00379 for (unsigned int i=0; i<d2phi.size(); i++) 00380 for (unsigned int p=0; p<d2phi[i].size(); p++) 00381 { 00382 d2phi[i][p].slice(0).slice(0) = d2phidx2[i][p] = 00383 d2phidxi2[i][p]*dxidx_map[p]*dxidx_map[p] + 00384 2*d2phidxideta[i][p]*dxidx_map[p]*detadx_map[p] + 00385 2*d2phidxidzeta[i][p]*dxidx_map[p]*dzetadx_map[p] + 00386 2*d2phidetadzeta[i][p]*detadx_map[p]*dzetadx_map[p] + 00387 d2phideta2[i][p]*detadx_map[p]*detadx_map[p] + 00388 d2phidzeta2[i][p]*dzetadx_map[p]*dzetadx_map[p]; 00389 00390 d2phi[i][p].slice(0).slice(1) = 00391 d2phi[i][p].slice(1).slice(0) = d2phidxdy[i][p] = 00392 d2phidxi2[i][p]*dxidx_map[p]*dxidy_map[p] + 00393 d2phidxideta[i][p]*dxidx_map[p]*detady_map[p] + 00394 d2phidxidzeta[i][p]*dxidx_map[p]*dzetady_map[p] + 00395 d2phideta2[i][p]*detadx_map[p]*detady_map[p] + 00396 d2phidxideta[i][p]*detadx_map[p]*dxidy_map[p] + 00397 d2phidetadzeta[i][p]*detadx_map[p]*dzetady_map[p] + 00398 d2phidzeta2[i][p]*dzetadx_map[p]*dzetady_map[p] + 00399 d2phidxidzeta[i][p]*dzetadx_map[p]*dxidy_map[p] + 00400 d2phidetadzeta[i][p]*dzetadx_map[p]*detady_map[p]; 00401 00402 d2phi[i][p].slice(0).slice(2) = 00403 d2phi[i][p].slice(2).slice(0) = d2phidy2[i][p] = 00404 d2phidxi2[i][p]*dxidx_map[p]*dxidz_map[p] + 00405 d2phidxideta[i][p]*dxidx_map[p]*detadz_map[p] + 00406 d2phidxidzeta[i][p]*dxidx_map[p]*dzetadz_map[p] + 00407 d2phideta2[i][p]*detadx_map[p]*detadz_map[p] + 00408 d2phidxideta[i][p]*detadx_map[p]*dxidz_map[p] + 00409 d2phidetadzeta[i][p]*detadx_map[p]*dzetadz_map[p] + 00410 d2phidzeta2[i][p]*dzetadx_map[p]*dzetadz_map[p] + 00411 d2phidxidzeta[i][p]*dzetadx_map[p]*dxidz_map[p] + 00412 d2phidetadzeta[i][p]*dzetadx_map[p]*detadz_map[p]; 00413 00414 d2phi[i][p].slice(1).slice(1) = d2phidxdz[i][p] = 00415 d2phidxi2[i][p]*dxidy_map[p]*dxidy_map[p] + 00416 2*d2phidxideta[i][p]*dxidy_map[p]*detady_map[p] + 00417 2*d2phidxidzeta[i][p]*dxidy_map[p]*dzetady_map[p] + 00418 2*d2phidetadzeta[i][p]*detady_map[p]*dzetady_map[p] + 00419 d2phideta2[i][p]*detady_map[p]*detady_map[p] + 00420 d2phidzeta2[i][p]*dzetady_map[p]*dzetady_map[p]; 00421 00422 d2phi[i][p].slice(1).slice(2) = 00423 d2phi[i][p].slice(2).slice(1) = d2phidydz[i][p] = 00424 d2phidxi2[i][p]*dxidy_map[p]*dxidz_map[p] + 00425 d2phidxideta[i][p]*dxidy_map[p]*detadz_map[p] + 00426 d2phidxidzeta[i][p]*dxidy_map[p]*dzetadz_map[p] + 00427 d2phideta2[i][p]*detady_map[p]*detadz_map[p] + 00428 d2phidxideta[i][p]*detady_map[p]*dxidz_map[p] + 00429 d2phidetadzeta[i][p]*detady_map[p]*dzetadz_map[p] + 00430 d2phidzeta2[i][p]*dzetady_map[p]*dzetadz_map[p] + 00431 d2phidxidzeta[i][p]*dzetady_map[p]*dxidz_map[p] + 00432 d2phidetadzeta[i][p]*dzetady_map[p]*detadz_map[p]; 00433 00434 d2phi[i][p].slice(2).slice(2) = d2phidz2[i][p] = 00435 d2phidxi2[i][p]*dxidz_map[p]*dxidz_map[p] + 00436 2*d2phidxideta[i][p]*dxidz_map[p]*detadz_map[p] + 00437 2*d2phidxidzeta[i][p]*dxidz_map[p]*dzetadz_map[p] + 00438 2*d2phidetadzeta[i][p]*detadz_map[p]*dzetadz_map[p] + 00439 d2phideta2[i][p]*detadz_map[p]*detadz_map[p] + 00440 d2phidzeta2[i][p]*dzetadz_map[p]*dzetadz_map[p]; 00441 } 00442 00443 break; 00444 } 00445 00446 default: 00447 libmesh_error(); 00448 00449 } // switch(dim) 00450 00451 return; 00452 } 00453 #endif //LIBMESH_ENABLE_SECOND_DERIVATIVES 00454 00455 template< > 00456 void H1FETransformation<Real>::map_curl( const unsigned int, 00457 const Elem* const, 00458 const std::vector<Point>&, 00459 const FEGenericBase<Real>&, 00460 std::vector<std::vector<Real> >& ) const 00461 { 00462 libMesh::err << "Computing the curl of a shape function only\n" 00463 << "makes sense for vector-valued elements." << std::endl; 00464 libmesh_error(); 00465 } 00466 00467 template< > 00468 void H1FETransformation<RealGradient>::map_curl( const unsigned int dim, 00469 const Elem* const, 00470 const std::vector<Point>&, 00471 const FEGenericBase<RealGradient>& fe, 00472 std::vector<std::vector<RealGradient> >& curl_phi ) const 00473 { 00474 switch(dim) 00475 { 00476 // The curl only make sense in 2D and 3D 00477 case 0: 00478 case 1: 00479 { 00480 libmesh_error(); 00481 } 00482 case 2: 00483 { 00484 const std::vector<std::vector<RealGradient> >& dphidxi = fe.get_dphidxi(); 00485 const std::vector<std::vector<RealGradient> >& dphideta = fe.get_dphideta(); 00486 00487 const std::vector<Real>& dxidx_map = fe.get_fe_map().get_dxidx(); 00488 const std::vector<Real>& dxidy_map = fe.get_fe_map().get_dxidy(); 00489 #if LIBMESH_DIM > 2 00490 const std::vector<Real>& dxidz_map = fe.get_fe_map().get_dxidz(); 00491 #endif 00492 00493 const std::vector<Real>& detadx_map = fe.get_fe_map().get_detadx(); 00494 const std::vector<Real>& detady_map = fe.get_fe_map().get_detady(); 00495 #if LIBMESH_DIM > 2 00496 const std::vector<Real>& detadz_map = fe.get_fe_map().get_detadz(); 00497 #endif 00498 00499 /* 00500 For 2D elements in 3D space: 00501 00502 curl = ( -dphi_y/dz, dphi_x/dz, dphi_y/dx - dphi_x/dy ) 00503 */ 00504 for (unsigned int i=0; i<curl_phi.size(); i++) 00505 for (unsigned int p=0; p<curl_phi[i].size(); p++) 00506 { 00507 00508 Real dphiy_dx = (dphidxi[i][p].slice(1))*dxidx_map[p] 00509 + (dphideta[i][p].slice(1))*detadx_map[p]; 00510 00511 Real dphix_dy = (dphidxi[i][p].slice(0))*dxidy_map[p] 00512 + (dphideta[i][p].slice(0))*detady_map[p]; 00513 00514 curl_phi[i][p].slice(2) = dphiy_dx - dphix_dy; 00515 00516 #if LIBMESH_DIM > 2 00517 Real dphiy_dz = (dphidxi[i][p].slice(1))*dxidz_map[p] 00518 + (dphideta[i][p].slice(1))*detadz_map[p]; 00519 00520 Real dphix_dz = (dphidxi[i][p].slice(0))*dxidz_map[p] 00521 + (dphideta[i][p].slice(0))*detadz_map[p]; 00522 00523 curl_phi[i][p].slice(0) = -dphiy_dz; 00524 curl_phi[i][p].slice(1) = dphix_dz; 00525 #endif 00526 } 00527 00528 break; 00529 } 00530 case 3: 00531 { 00532 const std::vector<std::vector<RealGradient> >& dphidxi = fe.get_dphidxi(); 00533 const std::vector<std::vector<RealGradient> >& dphideta = fe.get_dphideta(); 00534 const std::vector<std::vector<RealGradient> >& dphidzeta = fe.get_dphidzeta(); 00535 00536 const std::vector<Real>& dxidx_map = fe.get_fe_map().get_dxidx(); 00537 const std::vector<Real>& dxidy_map = fe.get_fe_map().get_dxidy(); 00538 const std::vector<Real>& dxidz_map = fe.get_fe_map().get_dxidz(); 00539 00540 const std::vector<Real>& detadx_map = fe.get_fe_map().get_detadx(); 00541 const std::vector<Real>& detady_map = fe.get_fe_map().get_detady(); 00542 const std::vector<Real>& detadz_map = fe.get_fe_map().get_detadz(); 00543 00544 const std::vector<Real>& dzetadx_map = fe.get_fe_map().get_dzetadx(); 00545 const std::vector<Real>& dzetady_map = fe.get_fe_map().get_dzetady(); 00546 const std::vector<Real>& dzetadz_map = fe.get_fe_map().get_dzetadz(); 00547 00548 /* 00549 In 3D: curl = ( dphi_z/dy - dphi_y/dz, dphi_x/dz - dphi_z/dx, dphi_y/dx - dphi_x/dy ) 00550 */ 00551 for (unsigned int i=0; i<curl_phi.size(); i++) 00552 for (unsigned int p=0; p<curl_phi[i].size(); p++) 00553 { 00554 Real dphiz_dy = (dphidxi[i][p].slice(2))*dxidy_map[p] 00555 + (dphideta[i][p].slice(2))*detady_map[p] 00556 + (dphidzeta[i][p].slice(2))*dzetady_map[p]; 00557 00558 Real dphiy_dz = (dphidxi[i][p].slice(1))*dxidz_map[p] 00559 + (dphideta[i][p].slice(1))*detadz_map[p] 00560 + (dphidzeta[i][p].slice(1))*dzetadz_map[p]; 00561 00562 Real dphix_dz = (dphidxi[i][p].slice(0))*dxidz_map[p] 00563 + (dphideta[i][p].slice(0))*detadz_map[p] 00564 + (dphidzeta[i][p].slice(0))*dzetadz_map[p]; 00565 00566 Real dphiz_dx = (dphidxi[i][p].slice(2))*dxidx_map[p] 00567 + (dphideta[i][p].slice(2))*detadx_map[p] 00568 + (dphidzeta[i][p].slice(2))*dzetadx_map[p]; 00569 00570 Real dphiy_dx = (dphidxi[i][p].slice(1))*dxidx_map[p] 00571 + (dphideta[i][p].slice(1))*detadx_map[p] 00572 + (dphidzeta[i][p].slice(1))*dzetadx_map[p]; 00573 00574 Real dphix_dy = (dphidxi[i][p].slice(0))*dxidy_map[p] 00575 + (dphideta[i][p].slice(0))*detady_map[p] 00576 + (dphidzeta[i][p].slice(0))*dzetady_map[p]; 00577 00578 curl_phi[i][p].slice(0) = dphiz_dy - dphiy_dz; 00579 00580 curl_phi[i][p].slice(1) = dphix_dz - dphiz_dx; 00581 00582 curl_phi[i][p].slice(2) = dphiy_dx - dphix_dy; 00583 } 00584 00585 break; 00586 } 00587 default: 00588 libmesh_error(); 00589 } 00590 00591 return; 00592 } 00593 00594 00595 template< > 00596 void H1FETransformation<Real>::map_div 00597 (const unsigned int, 00598 const Elem* const, 00599 const std::vector<Point>&, 00600 const FEGenericBase<Real>&, 00601 std::vector<std::vector<FEGenericBase<Real>::OutputDivergence> >& ) const 00602 { 00603 libMesh::err << "Computing the divergence of a shape function only\n" 00604 << "makes sense for vector-valued elements." << std::endl; 00605 libmesh_error(); 00606 } 00607 00608 00609 template<> 00610 void H1FETransformation<RealGradient>::map_div 00611 (const unsigned int dim, 00612 const Elem* const, 00613 const std::vector<Point>&, 00614 const FEGenericBase<RealGradient>& fe, 00615 std::vector<std::vector<FEGenericBase<RealGradient>::OutputDivergence> >& div_phi) const 00616 { 00617 switch(dim) 00618 { 00619 // The divergence only make sense in 2D and 3D 00620 case 0: 00621 case 1: 00622 { 00623 libmesh_error(); 00624 } 00625 case 2: 00626 { 00627 const std::vector<std::vector<RealGradient> >& dphidxi = fe.get_dphidxi(); 00628 const std::vector<std::vector<RealGradient> >& dphideta = fe.get_dphideta(); 00629 00630 const std::vector<Real>& dxidx_map = fe.get_fe_map().get_dxidx(); 00631 const std::vector<Real>& dxidy_map = fe.get_fe_map().get_dxidy(); 00632 00633 const std::vector<Real>& detadx_map = fe.get_fe_map().get_detadx(); 00634 const std::vector<Real>& detady_map = fe.get_fe_map().get_detady(); 00635 00636 /* 00637 In 2D: div = dphi_x/dx + dphi_y/dy 00638 */ 00639 for (unsigned int i=0; i<div_phi.size(); i++) 00640 for (unsigned int p=0; p<div_phi[i].size(); p++) 00641 { 00642 00643 Real dphix_dx = (dphidxi[i][p].slice(0))*dxidx_map[p] 00644 + (dphideta[i][p].slice(0))*detadx_map[p]; 00645 00646 Real dphiy_dy = (dphidxi[i][p].slice(1))*dxidy_map[p] 00647 + (dphideta[i][p].slice(1))*detady_map[p]; 00648 00649 div_phi[i][p] = dphix_dx + dphiy_dy; 00650 } 00651 break; 00652 } 00653 case 3: 00654 { 00655 const std::vector<std::vector<RealGradient> >& dphidxi = fe.get_dphidxi(); 00656 const std::vector<std::vector<RealGradient> >& dphideta = fe.get_dphideta(); 00657 const std::vector<std::vector<RealGradient> >& dphidzeta = fe.get_dphidzeta(); 00658 00659 const std::vector<Real>& dxidx_map = fe.get_fe_map().get_dxidx(); 00660 const std::vector<Real>& dxidy_map = fe.get_fe_map().get_dxidy(); 00661 const std::vector<Real>& dxidz_map = fe.get_fe_map().get_dxidz(); 00662 00663 const std::vector<Real>& detadx_map = fe.get_fe_map().get_detadx(); 00664 const std::vector<Real>& detady_map = fe.get_fe_map().get_detady(); 00665 const std::vector<Real>& detadz_map = fe.get_fe_map().get_detadz(); 00666 00667 const std::vector<Real>& dzetadx_map = fe.get_fe_map().get_dzetadx(); 00668 const std::vector<Real>& dzetady_map = fe.get_fe_map().get_dzetady(); 00669 const std::vector<Real>& dzetadz_map = fe.get_fe_map().get_dzetadz(); 00670 00671 /* 00672 In 3D: div = dphi_x/dx + dphi_y/dy + dphi_z/dz 00673 */ 00674 for (unsigned int i=0; i<div_phi.size(); i++) 00675 for (unsigned int p=0; p<div_phi[i].size(); p++) 00676 { 00677 Real dphix_dx = (dphidxi[i][p].slice(0))*dxidx_map[p] 00678 + (dphideta[i][p].slice(0))*detadx_map[p] 00679 + (dphidzeta[i][p].slice(0))*dzetadx_map[p]; 00680 00681 Real dphiy_dy = (dphidxi[i][p].slice(1))*dxidy_map[p] 00682 + (dphideta[i][p].slice(1))*detady_map[p] 00683 + (dphidzeta[i][p].slice(1))*dzetady_map[p]; 00684 00685 Real dphiz_dz = (dphidxi[i][p].slice(2))*dxidz_map[p] 00686 + (dphideta[i][p].slice(2))*detadz_map[p] 00687 + (dphidzeta[i][p].slice(2))*dzetadz_map[p]; 00688 00689 div_phi[i][p] = dphix_dx + dphiy_dy + dphiz_dz; 00690 } 00691 00692 break; 00693 } 00694 } // switch(dim) 00695 00696 return; 00697 } 00698 00699 // Explicit Instantations 00700 template class H1FETransformation<Real>; 00701 template class H1FETransformation<RealGradient>; 00702 00703 } //namespace libMesh
Site Created By: libMesh Developers
Last modified: February 05 2013 19:54:47 UTC
Hosted By: