h1_fe_transformation.C
Go to the documentation of this file.
00001 // The libMesh Finite Element Library. 00002 // Copyright (C) 2002-2014 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 00074 default: 00075 libmesh_error_msg("Invalid dim = " << dim); 00076 } 00077 } 00078 00079 00080 template< typename OutputShape > 00081 void H1FETransformation<OutputShape>::map_dphi( const unsigned int dim, 00082 const Elem* const, 00083 const std::vector<Point>&, 00084 const FEGenericBase<OutputShape>& fe, 00085 std::vector<std::vector<typename FEGenericBase<OutputShape>::OutputGradient> >& dphi, 00086 std::vector<std::vector<OutputShape> >& dphidx, 00087 std::vector<std::vector<OutputShape> >& dphidy, 00088 std::vector<std::vector<OutputShape> >& dphidz ) const 00089 { 00090 switch(dim) 00091 { 00092 case 0: // No derivatives in 0D 00093 { 00094 for (unsigned int i=0; i<dphi.size(); i++) 00095 for (unsigned int p=0; p<dphi[i].size(); p++) 00096 { 00097 dphi[i][p] = 0.; 00098 } 00099 break; 00100 } 00101 00102 case 1: 00103 { 00104 const std::vector<std::vector<OutputShape> >& dphidxi = fe.get_dphidxi(); 00105 00106 const std::vector<Real>& dxidx_map = fe.get_fe_map().get_dxidx(); 00107 #if LIBMESH_DIM>1 00108 const std::vector<Real>& dxidy_map = fe.get_fe_map().get_dxidy(); 00109 #endif 00110 #if LIBMESH_DIM>2 00111 const std::vector<Real>& dxidz_map = fe.get_fe_map().get_dxidz(); 00112 #endif 00113 00114 for (unsigned int i=0; i<dphi.size(); i++) 00115 for (unsigned int p=0; p<dphi[i].size(); p++) 00116 { 00117 // dphi/dx = (dphi/dxi)*(dxi/dx) 00118 dphi[i][p].slice(0) = dphidx[i][p] = dphidxi[i][p]*dxidx_map[p]; 00119 00120 #if LIBMESH_DIM>1 00121 dphi[i][p].slice(1) = dphidy[i][p] = dphidxi[i][p]*dxidy_map[p]; 00122 #endif 00123 #if LIBMESH_DIM>2 00124 dphi[i][p].slice(2) = dphidz[i][p] = dphidxi[i][p]*dxidz_map[p]; 00125 #endif 00126 } 00127 00128 break; 00129 } 00130 00131 case 2: 00132 { 00133 const std::vector<std::vector<OutputShape> >& dphidxi = fe.get_dphidxi(); 00134 const std::vector<std::vector<OutputShape> >& dphideta = fe.get_dphideta(); 00135 00136 const std::vector<Real>& dxidx_map = fe.get_fe_map().get_dxidx(); 00137 const std::vector<Real>& dxidy_map = fe.get_fe_map().get_dxidy(); 00138 #if LIBMESH_DIM > 2 00139 const std::vector<Real>& dxidz_map = fe.get_fe_map().get_dxidz(); 00140 #endif 00141 00142 const std::vector<Real>& detadx_map = fe.get_fe_map().get_detadx(); 00143 const std::vector<Real>& detady_map = fe.get_fe_map().get_detady(); 00144 #if LIBMESH_DIM > 2 00145 const std::vector<Real>& detadz_map = fe.get_fe_map().get_detadz(); 00146 #endif 00147 00148 for (unsigned int i=0; i<dphi.size(); i++) 00149 for (unsigned int p=0; p<dphi[i].size(); p++) 00150 { 00151 // dphi/dx = (dphi/dxi)*(dxi/dx) + (dphi/deta)*(deta/dx) 00152 dphi[i][p].slice(0) = dphidx[i][p] = (dphidxi[i][p]*dxidx_map[p] + 00153 dphideta[i][p]*detadx_map[p]); 00154 00155 // dphi/dy = (dphi/dxi)*(dxi/dy) + (dphi/deta)*(deta/dy) 00156 dphi[i][p].slice(1) = dphidy[i][p] = (dphidxi[i][p]*dxidy_map[p] + 00157 dphideta[i][p]*detady_map[p]); 00158 00159 #if LIBMESH_DIM > 2 00160 // dphi/dz = (dphi/dxi)*(dxi/dz) + (dphi/deta)*(deta/dz) 00161 dphi[i][p].slice(2) = dphidz[i][p] = (dphidxi[i][p]*dxidz_map[p] + 00162 dphideta[i][p]*detadz_map[p]); 00163 #endif 00164 } 00165 00166 break; 00167 } 00168 00169 case 3: 00170 { 00171 const std::vector<std::vector<OutputShape> >& dphidxi = fe.get_dphidxi(); 00172 const std::vector<std::vector<OutputShape> >& dphideta = fe.get_dphideta(); 00173 const std::vector<std::vector<OutputShape> >& dphidzeta = fe.get_dphidzeta(); 00174 00175 const std::vector<Real>& dxidx_map = fe.get_fe_map().get_dxidx(); 00176 const std::vector<Real>& dxidy_map = fe.get_fe_map().get_dxidy(); 00177 const std::vector<Real>& dxidz_map = fe.get_fe_map().get_dxidz(); 00178 00179 const std::vector<Real>& detadx_map = fe.get_fe_map().get_detadx(); 00180 const std::vector<Real>& detady_map = fe.get_fe_map().get_detady(); 00181 const std::vector<Real>& detadz_map = fe.get_fe_map().get_detadz(); 00182 00183 const std::vector<Real>& dzetadx_map = fe.get_fe_map().get_dzetadx(); 00184 const std::vector<Real>& dzetady_map = fe.get_fe_map().get_dzetady(); 00185 const std::vector<Real>& dzetadz_map = fe.get_fe_map().get_dzetadz(); 00186 00187 for (unsigned int i=0; i<dphi.size(); i++) 00188 for (unsigned int p=0; p<dphi[i].size(); p++) 00189 { 00190 // dphi/dx = (dphi/dxi)*(dxi/dx) + (dphi/deta)*(deta/dx) + (dphi/dzeta)*(dzeta/dx); 00191 dphi[i][p].slice(0) = dphidx[i][p] = (dphidxi[i][p]*dxidx_map[p] + 00192 dphideta[i][p]*detadx_map[p] + 00193 dphidzeta[i][p]*dzetadx_map[p]); 00194 00195 // dphi/dy = (dphi/dxi)*(dxi/dy) + (dphi/deta)*(deta/dy) + (dphi/dzeta)*(dzeta/dy); 00196 dphi[i][p].slice(1) = dphidy[i][p] = (dphidxi[i][p]*dxidy_map[p] + 00197 dphideta[i][p]*detady_map[p] + 00198 dphidzeta[i][p]*dzetady_map[p]); 00199 00200 // dphi/dz = (dphi/dxi)*(dxi/dz) + (dphi/deta)*(deta/dz) + (dphi/dzeta)*(dzeta/dz); 00201 dphi[i][p].slice(2) = dphidz[i][p] = (dphidxi[i][p]*dxidz_map[p] + 00202 dphideta[i][p]*detadz_map[p] + 00203 dphidzeta[i][p]*dzetadz_map[p]); 00204 } 00205 break; 00206 } 00207 00208 default: 00209 libmesh_error_msg("Invalid dim = " << dim); 00210 } // switch(dim) 00211 } 00212 00213 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES 00214 template< typename OutputShape > 00215 void H1FETransformation<OutputShape>::map_d2phi( const unsigned int dim, 00216 const Elem* const elem, 00217 const std::vector<Point>&, 00218 const FEGenericBase<OutputShape>& fe, 00219 std::vector<std::vector<typename FEGenericBase<OutputShape>::OutputTensor> >& d2phi, 00220 std::vector<std::vector<OutputShape> >& d2phidx2, 00221 std::vector<std::vector<OutputShape> >& d2phidxdy, 00222 std::vector<std::vector<OutputShape> >& d2phidxdz, 00223 std::vector<std::vector<OutputShape> >& d2phidy2, 00224 std::vector<std::vector<OutputShape> >& d2phidydz, 00225 std::vector<std::vector<OutputShape> >& d2phidz2 ) const 00226 { 00227 libmesh_do_once( 00228 if (!elem->has_affine_map()) 00229 { 00230 libMesh::err << "WARNING: Second derivatives are not currently " 00231 << "correctly calculated on non-affine elements!" 00232 << std::endl; 00233 } 00234 ); 00235 00236 00237 switch (dim) 00238 { 00239 case 0: // No derivatives in 0D 00240 { 00241 for (unsigned int i=0; i<d2phi.size(); i++) 00242 for (unsigned int p=0; p<d2phi[i].size(); p++) 00243 { 00244 d2phi[i][p] = 0.; 00245 } 00246 00247 break; 00248 } 00249 00250 case 1: 00251 { 00252 const std::vector<std::vector<OutputShape> >& d2phidxi2 = fe.get_d2phidxi2(); 00253 00254 const std::vector<Real>& dxidx_map = fe.get_fe_map().get_dxidx(); 00255 #if LIBMESH_DIM>1 00256 const std::vector<Real>& dxidy_map = fe.get_fe_map().get_dxidy(); 00257 #endif 00258 #if LIBMESH_DIM>2 00259 const std::vector<Real>& dxidz_map = fe.get_fe_map().get_dxidz(); 00260 #endif 00261 00262 for (unsigned int i=0; i<d2phi.size(); i++) 00263 for (unsigned int p=0; p<d2phi[i].size(); p++) 00264 { 00265 d2phi[i][p].slice(0).slice(0) = d2phidx2[i][p] = 00266 d2phidxi2[i][p]*dxidx_map[p]*dxidx_map[p]; 00267 #if LIBMESH_DIM>1 00268 d2phi[i][p].slice(0).slice(1) = 00269 d2phi[i][p].slice(1).slice(0) = d2phidxdy[i][p] = 00270 d2phidxi2[i][p]*dxidx_map[p]*dxidy_map[p]; 00271 00272 d2phi[i][p].slice(1).slice(1) = d2phidy2[i][p] = 00273 d2phidxi2[i][p]*dxidy_map[p]*dxidy_map[p]; 00274 #endif 00275 #if LIBMESH_DIM>2 00276 d2phi[i][p].slice(0).slice(2) = 00277 d2phi[i][p].slice(2).slice(0) = d2phidxdz[i][p] = 00278 d2phidxi2[i][p]*dxidx_map[p]*dxidz_map[p]; 00279 00280 d2phi[i][p].slice(1).slice(2) = 00281 d2phi[i][p].slice(2).slice(1) = d2phidydz[i][p] = 00282 d2phidxi2[i][p]*dxidy_map[p]*dxidz_map[p]; 00283 00284 d2phi[i][p].slice(2).slice(2) = d2phidz2[i][p] = 00285 d2phidxi2[i][p]*dxidz_map[p]*dxidz_map[p]; 00286 #endif 00287 } 00288 break; 00289 } 00290 00291 case 2: 00292 { 00293 const std::vector<std::vector<OutputShape> >& d2phidxi2 = fe.get_d2phidxi2(); 00294 const std::vector<std::vector<OutputShape> >& d2phidxideta = fe.get_d2phidxideta(); 00295 const std::vector<std::vector<OutputShape> >& d2phideta2 = fe.get_d2phideta2(); 00296 00297 const std::vector<Real>& dxidx_map = fe.get_fe_map().get_dxidx(); 00298 const std::vector<Real>& dxidy_map = fe.get_fe_map().get_dxidy(); 00299 #if LIBMESH_DIM > 2 00300 const std::vector<Real>& dxidz_map = fe.get_fe_map().get_dxidz(); 00301 #endif 00302 00303 const std::vector<Real>& detadx_map = fe.get_fe_map().get_detadx(); 00304 const std::vector<Real>& detady_map = fe.get_fe_map().get_detady(); 00305 #if LIBMESH_DIM > 2 00306 const std::vector<Real>& detadz_map = fe.get_fe_map().get_detadz(); 00307 #endif 00308 00309 for (unsigned int i=0; i<d2phi.size(); i++) 00310 for (unsigned int p=0; p<d2phi[i].size(); p++) 00311 { 00312 d2phi[i][p].slice(0).slice(0) = d2phidx2[i][p] = 00313 d2phidxi2[i][p]*dxidx_map[p]*dxidx_map[p] + 00314 2*d2phidxideta[i][p]*dxidx_map[p]*detadx_map[p] + 00315 d2phideta2[i][p]*detadx_map[p]*detadx_map[p]; 00316 00317 d2phi[i][p].slice(0).slice(1) = 00318 d2phi[i][p].slice(1).slice(0) = d2phidxdy[i][p] = 00319 d2phidxi2[i][p]*dxidx_map[p]*dxidy_map[p] + 00320 d2phidxideta[i][p]*dxidx_map[p]*detady_map[p] + 00321 d2phideta2[i][p]*detadx_map[p]*detady_map[p] + 00322 d2phidxideta[i][p]*detadx_map[p]*dxidy_map[p]; 00323 00324 d2phi[i][p].slice(1).slice(1) = d2phidy2[i][p] = 00325 d2phidxi2[i][p]*dxidy_map[p]*dxidy_map[p] + 00326 2*d2phidxideta[i][p]*dxidy_map[p]*detady_map[p] + 00327 d2phideta2[i][p]*detady_map[p]*detady_map[p]; 00328 00329 #if LIBMESH_DIM > 2 00330 d2phi[i][p].slice(0).slice(2) = 00331 d2phi[i][p].slice(2).slice(0) = d2phidxdz[i][p] = 00332 d2phidxi2[i][p]*dxidx_map[p]*dxidz_map[p] + 00333 d2phidxideta[i][p]*dxidx_map[p]*detadz_map[p] + 00334 d2phideta2[i][p]*detadx_map[p]*detadz_map[p] + 00335 d2phidxideta[i][p]*detadx_map[p]*dxidz_map[p]; 00336 00337 d2phi[i][p].slice(1).slice(2) = 00338 d2phi[i][p].slice(2).slice(1) = d2phidydz[i][p] = 00339 d2phidxi2[i][p]*dxidy_map[p]*dxidz_map[p] + 00340 d2phidxideta[i][p]*dxidy_map[p]*detadz_map[p] + 00341 d2phideta2[i][p]*detady_map[p]*detadz_map[p] + 00342 d2phidxideta[i][p]*detady_map[p]*dxidz_map[p]; 00343 00344 d2phi[i][p].slice(2).slice(2) = d2phidz2[i][p] = 00345 d2phidxi2[i][p]*dxidz_map[p]*dxidz_map[p] + 00346 2*d2phidxideta[i][p]*dxidz_map[p]*detadz_map[p] + 00347 d2phideta2[i][p]*detadz_map[p]*detadz_map[p]; 00348 #endif 00349 } 00350 00351 break; 00352 } 00353 00354 case 3: 00355 { 00356 const std::vector<std::vector<OutputShape> >& d2phidxi2 = fe.get_d2phidxi2(); 00357 const std::vector<std::vector<OutputShape> >& d2phidxideta = fe.get_d2phidxideta(); 00358 const std::vector<std::vector<OutputShape> >& d2phideta2 = fe.get_d2phideta2(); 00359 const std::vector<std::vector<OutputShape> >& d2phidxidzeta = fe.get_d2phidxidzeta(); 00360 const std::vector<std::vector<OutputShape> >& d2phidetadzeta = fe.get_d2phidetadzeta(); 00361 const std::vector<std::vector<OutputShape> >& d2phidzeta2 = fe.get_d2phidzeta2(); 00362 00363 const std::vector<Real>& dxidx_map = fe.get_fe_map().get_dxidx(); 00364 const std::vector<Real>& dxidy_map = fe.get_fe_map().get_dxidy(); 00365 const std::vector<Real>& dxidz_map = fe.get_fe_map().get_dxidz(); 00366 00367 const std::vector<Real>& detadx_map = fe.get_fe_map().get_detadx(); 00368 const std::vector<Real>& detady_map = fe.get_fe_map().get_detady(); 00369 const std::vector<Real>& detadz_map = fe.get_fe_map().get_detadz(); 00370 00371 const std::vector<Real>& dzetadx_map = fe.get_fe_map().get_dzetadx(); 00372 const std::vector<Real>& dzetady_map = fe.get_fe_map().get_dzetady(); 00373 const std::vector<Real>& dzetadz_map = fe.get_fe_map().get_dzetadz(); 00374 00375 for (unsigned int i=0; i<d2phi.size(); i++) 00376 for (unsigned int p=0; p<d2phi[i].size(); p++) 00377 { 00378 d2phi[i][p].slice(0).slice(0) = d2phidx2[i][p] = 00379 d2phidxi2[i][p]*dxidx_map[p]*dxidx_map[p] + 00380 2*d2phidxideta[i][p]*dxidx_map[p]*detadx_map[p] + 00381 2*d2phidxidzeta[i][p]*dxidx_map[p]*dzetadx_map[p] + 00382 2*d2phidetadzeta[i][p]*detadx_map[p]*dzetadx_map[p] + 00383 d2phideta2[i][p]*detadx_map[p]*detadx_map[p] + 00384 d2phidzeta2[i][p]*dzetadx_map[p]*dzetadx_map[p]; 00385 00386 d2phi[i][p].slice(0).slice(1) = 00387 d2phi[i][p].slice(1).slice(0) = d2phidxdy[i][p] = 00388 d2phidxi2[i][p]*dxidx_map[p]*dxidy_map[p] + 00389 d2phidxideta[i][p]*dxidx_map[p]*detady_map[p] + 00390 d2phidxidzeta[i][p]*dxidx_map[p]*dzetady_map[p] + 00391 d2phideta2[i][p]*detadx_map[p]*detady_map[p] + 00392 d2phidxideta[i][p]*detadx_map[p]*dxidy_map[p] + 00393 d2phidetadzeta[i][p]*detadx_map[p]*dzetady_map[p] + 00394 d2phidzeta2[i][p]*dzetadx_map[p]*dzetady_map[p] + 00395 d2phidxidzeta[i][p]*dzetadx_map[p]*dxidy_map[p] + 00396 d2phidetadzeta[i][p]*dzetadx_map[p]*detady_map[p]; 00397 00398 d2phi[i][p].slice(0).slice(2) = 00399 d2phi[i][p].slice(2).slice(0) = d2phidy2[i][p] = 00400 d2phidxi2[i][p]*dxidx_map[p]*dxidz_map[p] + 00401 d2phidxideta[i][p]*dxidx_map[p]*detadz_map[p] + 00402 d2phidxidzeta[i][p]*dxidx_map[p]*dzetadz_map[p] + 00403 d2phideta2[i][p]*detadx_map[p]*detadz_map[p] + 00404 d2phidxideta[i][p]*detadx_map[p]*dxidz_map[p] + 00405 d2phidetadzeta[i][p]*detadx_map[p]*dzetadz_map[p] + 00406 d2phidzeta2[i][p]*dzetadx_map[p]*dzetadz_map[p] + 00407 d2phidxidzeta[i][p]*dzetadx_map[p]*dxidz_map[p] + 00408 d2phidetadzeta[i][p]*dzetadx_map[p]*detadz_map[p]; 00409 00410 d2phi[i][p].slice(1).slice(1) = d2phidxdz[i][p] = 00411 d2phidxi2[i][p]*dxidy_map[p]*dxidy_map[p] + 00412 2*d2phidxideta[i][p]*dxidy_map[p]*detady_map[p] + 00413 2*d2phidxidzeta[i][p]*dxidy_map[p]*dzetady_map[p] + 00414 2*d2phidetadzeta[i][p]*detady_map[p]*dzetady_map[p] + 00415 d2phideta2[i][p]*detady_map[p]*detady_map[p] + 00416 d2phidzeta2[i][p]*dzetady_map[p]*dzetady_map[p]; 00417 00418 d2phi[i][p].slice(1).slice(2) = 00419 d2phi[i][p].slice(2).slice(1) = d2phidydz[i][p] = 00420 d2phidxi2[i][p]*dxidy_map[p]*dxidz_map[p] + 00421 d2phidxideta[i][p]*dxidy_map[p]*detadz_map[p] + 00422 d2phidxidzeta[i][p]*dxidy_map[p]*dzetadz_map[p] + 00423 d2phideta2[i][p]*detady_map[p]*detadz_map[p] + 00424 d2phidxideta[i][p]*detady_map[p]*dxidz_map[p] + 00425 d2phidetadzeta[i][p]*detady_map[p]*dzetadz_map[p] + 00426 d2phidzeta2[i][p]*dzetady_map[p]*dzetadz_map[p] + 00427 d2phidxidzeta[i][p]*dzetady_map[p]*dxidz_map[p] + 00428 d2phidetadzeta[i][p]*dzetady_map[p]*detadz_map[p]; 00429 00430 d2phi[i][p].slice(2).slice(2) = d2phidz2[i][p] = 00431 d2phidxi2[i][p]*dxidz_map[p]*dxidz_map[p] + 00432 2*d2phidxideta[i][p]*dxidz_map[p]*detadz_map[p] + 00433 2*d2phidxidzeta[i][p]*dxidz_map[p]*dzetadz_map[p] + 00434 2*d2phidetadzeta[i][p]*detadz_map[p]*dzetadz_map[p] + 00435 d2phideta2[i][p]*detadz_map[p]*detadz_map[p] + 00436 d2phidzeta2[i][p]*dzetadz_map[p]*dzetadz_map[p]; 00437 } 00438 00439 break; 00440 } 00441 00442 default: 00443 libmesh_error_msg("Invalid dim = " << dim); 00444 } // switch(dim) 00445 } 00446 #endif //LIBMESH_ENABLE_SECOND_DERIVATIVES 00447 00448 template< > 00449 void H1FETransformation<Real>::map_curl( const unsigned int, 00450 const Elem* const, 00451 const std::vector<Point>&, 00452 const FEGenericBase<Real>&, 00453 std::vector<std::vector<Real> >& ) const 00454 { 00455 libmesh_error_msg("Computing the curl of a shape function only \nmakes sense for vector-valued elements."); 00456 } 00457 00458 template< > 00459 void H1FETransformation<RealGradient>::map_curl( const unsigned int dim, 00460 const Elem* const, 00461 const std::vector<Point>&, 00462 const FEGenericBase<RealGradient>& fe, 00463 std::vector<std::vector<RealGradient> >& curl_phi ) const 00464 { 00465 switch (dim) 00466 { 00467 case 0: 00468 case 1: 00469 libmesh_error_msg("The curl only make sense in 2D and 3D"); 00470 00471 case 2: 00472 { 00473 const std::vector<std::vector<RealGradient> >& dphidxi = fe.get_dphidxi(); 00474 const std::vector<std::vector<RealGradient> >& dphideta = fe.get_dphideta(); 00475 00476 const std::vector<Real>& dxidx_map = fe.get_fe_map().get_dxidx(); 00477 const std::vector<Real>& dxidy_map = fe.get_fe_map().get_dxidy(); 00478 #if LIBMESH_DIM > 2 00479 const std::vector<Real>& dxidz_map = fe.get_fe_map().get_dxidz(); 00480 #endif 00481 00482 const std::vector<Real>& detadx_map = fe.get_fe_map().get_detadx(); 00483 const std::vector<Real>& detady_map = fe.get_fe_map().get_detady(); 00484 #if LIBMESH_DIM > 2 00485 const std::vector<Real>& detadz_map = fe.get_fe_map().get_detadz(); 00486 #endif 00487 00488 /* 00489 For 2D elements in 3D space: 00490 00491 curl = ( -dphi_y/dz, dphi_x/dz, dphi_y/dx - dphi_x/dy ) 00492 */ 00493 for (unsigned int i=0; i<curl_phi.size(); i++) 00494 for (unsigned int p=0; p<curl_phi[i].size(); p++) 00495 { 00496 00497 Real dphiy_dx = (dphidxi[i][p].slice(1))*dxidx_map[p] 00498 + (dphideta[i][p].slice(1))*detadx_map[p]; 00499 00500 Real dphix_dy = (dphidxi[i][p].slice(0))*dxidy_map[p] 00501 + (dphideta[i][p].slice(0))*detady_map[p]; 00502 00503 curl_phi[i][p].slice(2) = dphiy_dx - dphix_dy; 00504 00505 #if LIBMESH_DIM > 2 00506 Real dphiy_dz = (dphidxi[i][p].slice(1))*dxidz_map[p] 00507 + (dphideta[i][p].slice(1))*detadz_map[p]; 00508 00509 Real dphix_dz = (dphidxi[i][p].slice(0))*dxidz_map[p] 00510 + (dphideta[i][p].slice(0))*detadz_map[p]; 00511 00512 curl_phi[i][p].slice(0) = -dphiy_dz; 00513 curl_phi[i][p].slice(1) = dphix_dz; 00514 #endif 00515 } 00516 00517 break; 00518 } 00519 case 3: 00520 { 00521 const std::vector<std::vector<RealGradient> >& dphidxi = fe.get_dphidxi(); 00522 const std::vector<std::vector<RealGradient> >& dphideta = fe.get_dphideta(); 00523 const std::vector<std::vector<RealGradient> >& dphidzeta = fe.get_dphidzeta(); 00524 00525 const std::vector<Real>& dxidx_map = fe.get_fe_map().get_dxidx(); 00526 const std::vector<Real>& dxidy_map = fe.get_fe_map().get_dxidy(); 00527 const std::vector<Real>& dxidz_map = fe.get_fe_map().get_dxidz(); 00528 00529 const std::vector<Real>& detadx_map = fe.get_fe_map().get_detadx(); 00530 const std::vector<Real>& detady_map = fe.get_fe_map().get_detady(); 00531 const std::vector<Real>& detadz_map = fe.get_fe_map().get_detadz(); 00532 00533 const std::vector<Real>& dzetadx_map = fe.get_fe_map().get_dzetadx(); 00534 const std::vector<Real>& dzetady_map = fe.get_fe_map().get_dzetady(); 00535 const std::vector<Real>& dzetadz_map = fe.get_fe_map().get_dzetadz(); 00536 00537 /* 00538 In 3D: curl = ( dphi_z/dy - dphi_y/dz, dphi_x/dz - dphi_z/dx, dphi_y/dx - dphi_x/dy ) 00539 */ 00540 for (unsigned int i=0; i<curl_phi.size(); i++) 00541 for (unsigned int p=0; p<curl_phi[i].size(); p++) 00542 { 00543 Real dphiz_dy = (dphidxi[i][p].slice(2))*dxidy_map[p] 00544 + (dphideta[i][p].slice(2))*detady_map[p] 00545 + (dphidzeta[i][p].slice(2))*dzetady_map[p]; 00546 00547 Real dphiy_dz = (dphidxi[i][p].slice(1))*dxidz_map[p] 00548 + (dphideta[i][p].slice(1))*detadz_map[p] 00549 + (dphidzeta[i][p].slice(1))*dzetadz_map[p]; 00550 00551 Real dphix_dz = (dphidxi[i][p].slice(0))*dxidz_map[p] 00552 + (dphideta[i][p].slice(0))*detadz_map[p] 00553 + (dphidzeta[i][p].slice(0))*dzetadz_map[p]; 00554 00555 Real dphiz_dx = (dphidxi[i][p].slice(2))*dxidx_map[p] 00556 + (dphideta[i][p].slice(2))*detadx_map[p] 00557 + (dphidzeta[i][p].slice(2))*dzetadx_map[p]; 00558 00559 Real dphiy_dx = (dphidxi[i][p].slice(1))*dxidx_map[p] 00560 + (dphideta[i][p].slice(1))*detadx_map[p] 00561 + (dphidzeta[i][p].slice(1))*dzetadx_map[p]; 00562 00563 Real dphix_dy = (dphidxi[i][p].slice(0))*dxidy_map[p] 00564 + (dphideta[i][p].slice(0))*detady_map[p] 00565 + (dphidzeta[i][p].slice(0))*dzetady_map[p]; 00566 00567 curl_phi[i][p].slice(0) = dphiz_dy - dphiy_dz; 00568 00569 curl_phi[i][p].slice(1) = dphix_dz - dphiz_dx; 00570 00571 curl_phi[i][p].slice(2) = dphiy_dx - dphix_dy; 00572 } 00573 00574 break; 00575 } 00576 default: 00577 libmesh_error_msg("Invalid dim = " << dim); 00578 } 00579 } 00580 00581 00582 template< > 00583 void H1FETransformation<Real>::map_div 00584 (const unsigned int, 00585 const Elem* const, 00586 const std::vector<Point>&, 00587 const FEGenericBase<Real>&, 00588 std::vector<std::vector<FEGenericBase<Real>::OutputDivergence> >& ) const 00589 { 00590 libmesh_error_msg("Computing the divergence of a shape function only \nmakes sense for vector-valued elements."); 00591 } 00592 00593 00594 template<> 00595 void H1FETransformation<RealGradient>::map_div 00596 (const unsigned int dim, 00597 const Elem* const, 00598 const std::vector<Point>&, 00599 const FEGenericBase<RealGradient>& fe, 00600 std::vector<std::vector<FEGenericBase<RealGradient>::OutputDivergence> >& div_phi) const 00601 { 00602 switch (dim) 00603 { 00604 case 0: 00605 case 1: 00606 libmesh_error_msg("The divergence only make sense in 2D and 3D"); 00607 00608 case 2: 00609 { 00610 const std::vector<std::vector<RealGradient> >& dphidxi = fe.get_dphidxi(); 00611 const std::vector<std::vector<RealGradient> >& dphideta = fe.get_dphideta(); 00612 00613 const std::vector<Real>& dxidx_map = fe.get_fe_map().get_dxidx(); 00614 const std::vector<Real>& dxidy_map = fe.get_fe_map().get_dxidy(); 00615 00616 const std::vector<Real>& detadx_map = fe.get_fe_map().get_detadx(); 00617 const std::vector<Real>& detady_map = fe.get_fe_map().get_detady(); 00618 00619 /* 00620 In 2D: div = dphi_x/dx + dphi_y/dy 00621 */ 00622 for (unsigned int i=0; i<div_phi.size(); i++) 00623 for (unsigned int p=0; p<div_phi[i].size(); p++) 00624 { 00625 00626 Real dphix_dx = (dphidxi[i][p].slice(0))*dxidx_map[p] 00627 + (dphideta[i][p].slice(0))*detadx_map[p]; 00628 00629 Real dphiy_dy = (dphidxi[i][p].slice(1))*dxidy_map[p] 00630 + (dphideta[i][p].slice(1))*detady_map[p]; 00631 00632 div_phi[i][p] = dphix_dx + dphiy_dy; 00633 } 00634 break; 00635 } 00636 case 3: 00637 { 00638 const std::vector<std::vector<RealGradient> >& dphidxi = fe.get_dphidxi(); 00639 const std::vector<std::vector<RealGradient> >& dphideta = fe.get_dphideta(); 00640 const std::vector<std::vector<RealGradient> >& dphidzeta = fe.get_dphidzeta(); 00641 00642 const std::vector<Real>& dxidx_map = fe.get_fe_map().get_dxidx(); 00643 const std::vector<Real>& dxidy_map = fe.get_fe_map().get_dxidy(); 00644 const std::vector<Real>& dxidz_map = fe.get_fe_map().get_dxidz(); 00645 00646 const std::vector<Real>& detadx_map = fe.get_fe_map().get_detadx(); 00647 const std::vector<Real>& detady_map = fe.get_fe_map().get_detady(); 00648 const std::vector<Real>& detadz_map = fe.get_fe_map().get_detadz(); 00649 00650 const std::vector<Real>& dzetadx_map = fe.get_fe_map().get_dzetadx(); 00651 const std::vector<Real>& dzetady_map = fe.get_fe_map().get_dzetady(); 00652 const std::vector<Real>& dzetadz_map = fe.get_fe_map().get_dzetadz(); 00653 00654 /* 00655 In 3D: div = dphi_x/dx + dphi_y/dy + dphi_z/dz 00656 */ 00657 for (unsigned int i=0; i<div_phi.size(); i++) 00658 for (unsigned int p=0; p<div_phi[i].size(); p++) 00659 { 00660 Real dphix_dx = (dphidxi[i][p].slice(0))*dxidx_map[p] 00661 + (dphideta[i][p].slice(0))*detadx_map[p] 00662 + (dphidzeta[i][p].slice(0))*dzetadx_map[p]; 00663 00664 Real dphiy_dy = (dphidxi[i][p].slice(1))*dxidy_map[p] 00665 + (dphideta[i][p].slice(1))*detady_map[p] 00666 + (dphidzeta[i][p].slice(1))*dzetady_map[p]; 00667 00668 Real dphiz_dz = (dphidxi[i][p].slice(2))*dxidz_map[p] 00669 + (dphideta[i][p].slice(2))*detadz_map[p] 00670 + (dphidzeta[i][p].slice(2))*dzetadz_map[p]; 00671 00672 div_phi[i][p] = dphix_dx + dphiy_dy + dphiz_dz; 00673 } 00674 00675 break; 00676 } 00677 00678 default: \ 00679 libmesh_error_msg("Invalid dim = " << dim); \ 00680 } // switch(dim) 00681 00682 return; 00683 } 00684 00685 // Explicit Instantations 00686 template class H1FETransformation<Real>; 00687 template class H1FETransformation<RealGradient>; 00688 00689 } //namespace libMesh
generated by
