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