hcurl_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/hcurl_fe_transformation.h" 00019 #include "libmesh/fe_interface.h" 00020 00021 namespace libMesh 00022 { 00023 00024 template< typename OutputShape > 00025 void HCurlFETransformation<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 // These element transformations only make sense in 2D and 3D 00034 case 0: 00035 case 1: 00036 { 00037 libmesh_error(); 00038 } 00039 00040 case 2: 00041 { 00042 const std::vector<Real>& dxidx_map = fe.get_fe_map().get_dxidx(); 00043 const std::vector<Real>& dxidy_map = fe.get_fe_map().get_dxidy(); 00044 00045 const std::vector<Real>& detadx_map = fe.get_fe_map().get_detadx(); 00046 const std::vector<Real>& detady_map = fe.get_fe_map().get_detady(); 00047 00048 // FIXME: Need to update for 2D elements in 3D space 00049 /* phi = (dx/dxi)^-T * \hat{phi} 00050 In 2D: 00051 (dx/dxi)^{-1} = [ dxi/dx dxi/dy 00052 deta/dx deta/dy ] 00053 00054 so: dxi/dx^{-T} * \hat{phi} = [ dxi/dx deta/dx [ \hat{phi}_xi 00055 dxi/dy deta/dy ] \hat{phi}_eta ] 00056 00057 or in indicial notation: phi_j = xi_{i,j}*\hat{phi}_i */ 00058 00059 for (unsigned int i=0; i<phi.size(); i++) 00060 for (unsigned int p=0; p<phi[i].size(); p++) 00061 { 00062 // Need to temporarily cache reference shape functions 00063 // TODO: PB: Might be worth trying to build phi_ref separately to see 00064 // if we can get vectorization 00065 OutputShape phi_ref; 00066 FEInterface::shape<OutputShape>(2, fe.get_fe_type(), elem, i, qp[p], phi_ref); 00067 00068 phi[i][p](0) = dxidx_map[p]*phi_ref.slice(0) + detadx_map[p]*phi_ref.slice(1); 00069 00070 phi[i][p](1) = dxidy_map[p]*phi_ref.slice(0) + detady_map[p]*phi_ref.slice(1); 00071 } 00072 00073 break; 00074 } 00075 00076 case 3: 00077 { 00078 const std::vector<Real>& dxidx_map = fe.get_fe_map().get_dxidx(); 00079 const std::vector<Real>& dxidy_map = fe.get_fe_map().get_dxidy(); 00080 const std::vector<Real>& dxidz_map = fe.get_fe_map().get_dxidz(); 00081 00082 const std::vector<Real>& detadx_map = fe.get_fe_map().get_detadx(); 00083 const std::vector<Real>& detady_map = fe.get_fe_map().get_detady(); 00084 const std::vector<Real>& detadz_map = fe.get_fe_map().get_detadz(); 00085 00086 const std::vector<Real>& dzetadx_map = fe.get_fe_map().get_dzetadx(); 00087 const std::vector<Real>& dzetady_map = fe.get_fe_map().get_dzetady(); 00088 const std::vector<Real>& dzetadz_map = fe.get_fe_map().get_dzetadz(); 00089 00090 /* phi = (dx/dxi)^-T * \hat{phi} 00091 In 3D: 00092 dx/dxi^-1 = [ dxi/dx dxi/dy dxi/dz 00093 deta/dx deta/dy deta/dz 00094 dzeta/dx dzeta/dy dzeta/dz] 00095 00096 so: dxi/dx^-T * \hat{phi} = [ dxi/dx deta/dx dzeta/dx [ \hat{phi}_xi 00097 dxi/dy deta/dy dzeta/dy \hat{phi}_eta 00098 dxi/dz deta/dz dzeta/dz ] \hat{phi}_zeta ] 00099 00100 or in indicial notation: phi_j = xi_{i,j}*\hat{phi}_i */ 00101 00102 for (unsigned int i=0; i<phi.size(); i++) 00103 for (unsigned int p=0; p<phi[i].size(); p++) 00104 { 00105 // Need to temporarily cache reference shape functions 00106 // TODO: PB: Might be worth trying to build phi_ref separately to see 00107 // if we can get vectorization 00108 OutputShape phi_ref; 00109 FEInterface::shape<OutputShape>(3, fe.get_fe_type(), elem, i, qp[p], phi_ref); 00110 00111 phi[i][p].slice(0) = dxidx_map[p]*phi_ref.slice(0) + detadx_map[p]*phi_ref.slice(1) 00112 + dzetadx_map[p]*phi_ref.slice(2); 00113 00114 phi[i][p].slice(1) = dxidy_map[p]*phi_ref.slice(0) + detady_map[p]*phi_ref.slice(1) 00115 + dzetady_map[p]*phi_ref.slice(2); 00116 00117 phi[i][p].slice(2) = dxidz_map[p]*phi_ref.slice(0) + detadz_map[p]*phi_ref.slice(1) 00118 + dzetadz_map[p]*phi_ref.slice(2); 00119 } 00120 break; 00121 } 00122 00123 default: 00124 libmesh_error(); 00125 00126 } // switch(dim) 00127 00128 return; 00129 } 00130 00131 template< typename OutputShape > 00132 void HCurlFETransformation<OutputShape>::map_curl( const unsigned int dim, 00133 const Elem* const, 00134 const std::vector<Point>&, 00135 const FEGenericBase<OutputShape>& fe, 00136 std::vector<std::vector<OutputShape> >& curl_phi ) const 00137 { 00138 switch(dim) 00139 { 00140 // These element transformations only make sense in 2D and 3D 00141 case 0: 00142 case 1: 00143 { 00144 libmesh_error(); 00145 } 00146 00147 case 2: 00148 { 00149 const std::vector<std::vector<OutputShape> >& dphi_dxi = fe.get_dphidxi(); 00150 const std::vector<std::vector<OutputShape> >& dphi_deta = fe.get_dphideta(); 00151 00152 const std::vector<Real>& J = fe.get_fe_map().get_jacobian(); 00153 00154 // FIXME: I don't think this is valid for 2D elements in 3D space 00155 /* In 2D: curl(phi) = J^{-1} * curl(\hat{phi}) */ 00156 for (unsigned int i=0; i<curl_phi.size(); i++) 00157 for (unsigned int p=0; p<curl_phi[i].size(); p++) 00158 { 00159 curl_phi[i][p].slice(0) = curl_phi[i][p].slice(1) = 0.0; 00160 00161 curl_phi[i][p].slice(2) = ( dphi_dxi[i][p].slice(1) - dphi_deta[i][p].slice(0) )/J[p]; 00162 } 00163 00164 break; 00165 } 00166 00167 case 3: 00168 { 00169 const std::vector<std::vector<OutputShape> >& dphi_dxi = fe.get_dphidxi(); 00170 const std::vector<std::vector<OutputShape> >& dphi_deta = fe.get_dphideta(); 00171 const std::vector<std::vector<OutputShape> >& dphi_dzeta = fe.get_dphidzeta(); 00172 00173 const std::vector<RealGradient>& dxyz_dxi = fe.get_fe_map().get_dxyzdxi(); 00174 const std::vector<RealGradient>& dxyz_deta = fe.get_fe_map().get_dxyzdeta(); 00175 const std::vector<RealGradient>& dxyz_dzeta = fe.get_fe_map().get_dxyzdzeta(); 00176 00177 const std::vector<Real>& J = fe.get_fe_map().get_jacobian(); 00178 00179 for (unsigned int i=0; i<curl_phi.size(); i++) 00180 for (unsigned int p=0; p<curl_phi[i].size(); p++) 00181 { 00182 Real dx_dxi = dxyz_dxi[p](0); 00183 Real dx_deta = dxyz_deta[p](0); 00184 Real dx_dzeta = dxyz_dzeta[p](0); 00185 00186 Real dy_dxi = dxyz_dxi[p](1); 00187 Real dy_deta = dxyz_deta[p](1); 00188 Real dy_dzeta = dxyz_dzeta[p](1); 00189 00190 Real dz_dxi = dxyz_dxi[p](2); 00191 Real dz_deta = dxyz_deta[p](2); 00192 Real dz_dzeta = dxyz_dzeta[p](2); 00193 00194 const Real inv_jac = 1.0/J[p]; 00195 00196 /* In 3D: curl(phi) = J^{-1} dx/dxi * curl(\hat{phi}) 00197 00198 dx/dxi = [ dx/dxi dx/deta dx/dzeta 00199 dy/dxi dy/deta dy/dzeta 00200 dz/dxi dz/deta dz/dzeta ] 00201 00202 curl(u) = [ du_z/deta - du_y/dzeta 00203 du_x/dzeta - du_z/dxi 00204 du_y/dxi - du_x/deta ] 00205 */ 00206 curl_phi[i][p].slice(0) = inv_jac*( dx_dxi*( dphi_deta[i][p].slice(2) - 00207 dphi_dzeta[i][p].slice(1) ) + 00208 dx_deta*( dphi_dzeta[i][p].slice(0) - 00209 dphi_dxi[i][p].slice(2) ) + 00210 dx_dzeta*( dphi_dxi[i][p].slice(1) - 00211 dphi_deta[i][p].slice(0) ) ); 00212 00213 curl_phi[i][p].slice(1) = inv_jac*( dy_dxi*( dphi_deta[i][p].slice(2) - 00214 dphi_dzeta[i][p].slice(1) ) + 00215 dy_deta*( dphi_dzeta[i][p].slice(0)- 00216 dphi_dxi[i][p].slice(2) ) + 00217 dy_dzeta*( dphi_dxi[i][p].slice(1) - 00218 dphi_deta[i][p].slice(0) ) ); 00219 00220 curl_phi[i][p].slice(2) = inv_jac*( dz_dxi*( dphi_deta[i][p].slice(2) - 00221 dphi_dzeta[i][p].slice(1) ) + 00222 dz_deta*( dphi_dzeta[i][p].slice(0) - 00223 dphi_dxi[i][p].slice(2) ) + 00224 dz_dzeta*( dphi_dxi[i][p].slice(1) - 00225 dphi_deta[i][p].slice(0) ) ); 00226 } 00227 00228 break; 00229 } 00230 00231 default: 00232 libmesh_error(); 00233 00234 } // switch(dim) 00235 00236 return; 00237 } 00238 00239 template class HCurlFETransformation<RealGradient>; 00240 00241 template<> 00242 void HCurlFETransformation<Real>::map_phi( const unsigned int, 00243 const Elem* const, 00244 const std::vector<Point>&, 00245 const FEGenericBase<Real>&, 00246 std::vector<std::vector<Real> >& ) const 00247 { 00248 libMesh::err << "HCurl transformations only make sense for vector-valued elements." 00249 << std::endl; 00250 libmesh_error(); 00251 } 00252 00253 template<> 00254 void HCurlFETransformation<Real>::map_curl( const unsigned int, 00255 const Elem* const, 00256 const std::vector<Point>&, 00257 const FEGenericBase<Real>&, 00258 std::vector<std::vector<Real> >& ) const 00259 { 00260 libMesh::err << "HCurl transformations only make sense for vector-valued elements." 00261 << std::endl; 00262 libmesh_error(); 00263 } 00264 00265 00266 } // namespace libMesh
Site Created By: libMesh Developers
Last modified: February 05 2013 19:54:47 UTC
Hosted By: