hcurl_fe_transformation.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2014 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
3 
4 // This library is free software; you can redistribute it and/or
5 // modify it under the terms of the GNU Lesser General Public
6 // License as published by the Free Software Foundation; either
7 // version 2.1 of the License, or (at your option) any later version.
8 
9 // This library is distributed in the hope that it will be useful,
10 // but WITHOUT ANY WARRANTY; without even the implied warranty of
11 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12 // Lesser General Public License for more details.
13 
14 // You should have received a copy of the GNU Lesser General Public
15 // License along with this library; if not, write to the Free Software
16 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 
19 #include "libmesh/fe_interface.h"
20 
21 namespace libMesh
22 {
23 
24  template< typename OutputShape >
26  const Elem* const elem,
27  const std::vector<Point>& qp,
29  std::vector<std::vector<OutputShape> >& phi ) const
30  {
31  switch(dim)
32  {
33  // These element transformations only make sense in 2D and 3D
34  case 0:
35  case 1:
36  {
37  libmesh_error();
38  }
39 
40  case 2:
41  {
42  const std::vector<Real>& dxidx_map = fe.get_fe_map().get_dxidx();
43  const std::vector<Real>& dxidy_map = fe.get_fe_map().get_dxidy();
44 
45  const std::vector<Real>& detadx_map = fe.get_fe_map().get_detadx();
46  const std::vector<Real>& detady_map = fe.get_fe_map().get_detady();
47 
48  // FIXME: Need to update for 2D elements in 3D space
49  /* phi = (dx/dxi)^-T * \hat{phi}
50  In 2D:
51  (dx/dxi)^{-1} = [ dxi/dx dxi/dy
52  deta/dx deta/dy ]
53 
54  so: dxi/dx^{-T} * \hat{phi} = [ dxi/dx deta/dx [ \hat{phi}_xi
55  dxi/dy deta/dy ] \hat{phi}_eta ]
56 
57  or in indicial notation: phi_j = xi_{i,j}*\hat{phi}_i */
58 
59  for (unsigned int i=0; i<phi.size(); i++)
60  for (unsigned int p=0; p<phi[i].size(); p++)
61  {
62  // Need to temporarily cache reference shape functions
63  // TODO: PB: Might be worth trying to build phi_ref separately to see
64  // if we can get vectorization
65  OutputShape phi_ref;
66  FEInterface::shape<OutputShape>(2, fe.get_fe_type(), elem, i, qp[p], phi_ref);
67 
68  phi[i][p](0) = dxidx_map[p]*phi_ref.slice(0) + detadx_map[p]*phi_ref.slice(1);
69 
70  phi[i][p](1) = dxidy_map[p]*phi_ref.slice(0) + detady_map[p]*phi_ref.slice(1);
71  }
72 
73  break;
74  }
75 
76  case 3:
77  {
78  const std::vector<Real>& dxidx_map = fe.get_fe_map().get_dxidx();
79  const std::vector<Real>& dxidy_map = fe.get_fe_map().get_dxidy();
80  const std::vector<Real>& dxidz_map = fe.get_fe_map().get_dxidz();
81 
82  const std::vector<Real>& detadx_map = fe.get_fe_map().get_detadx();
83  const std::vector<Real>& detady_map = fe.get_fe_map().get_detady();
84  const std::vector<Real>& detadz_map = fe.get_fe_map().get_detadz();
85 
86  const std::vector<Real>& dzetadx_map = fe.get_fe_map().get_dzetadx();
87  const std::vector<Real>& dzetady_map = fe.get_fe_map().get_dzetady();
88  const std::vector<Real>& dzetadz_map = fe.get_fe_map().get_dzetadz();
89 
90  /* phi = (dx/dxi)^-T * \hat{phi}
91  In 3D:
92  dx/dxi^-1 = [ dxi/dx dxi/dy dxi/dz
93  deta/dx deta/dy deta/dz
94  dzeta/dx dzeta/dy dzeta/dz]
95 
96  so: dxi/dx^-T * \hat{phi} = [ dxi/dx deta/dx dzeta/dx [ \hat{phi}_xi
97  dxi/dy deta/dy dzeta/dy \hat{phi}_eta
98  dxi/dz deta/dz dzeta/dz ] \hat{phi}_zeta ]
99 
100  or in indicial notation: phi_j = xi_{i,j}*\hat{phi}_i */
101 
102  for (unsigned int i=0; i<phi.size(); i++)
103  for (unsigned int p=0; p<phi[i].size(); p++)
104  {
105  // Need to temporarily cache reference shape functions
106  // TODO: PB: Might be worth trying to build phi_ref separately to see
107  // if we can get vectorization
108  OutputShape phi_ref;
109  FEInterface::shape<OutputShape>(3, fe.get_fe_type(), elem, i, qp[p], phi_ref);
110 
111  phi[i][p].slice(0) = dxidx_map[p]*phi_ref.slice(0) + detadx_map[p]*phi_ref.slice(1)
112  + dzetadx_map[p]*phi_ref.slice(2);
113 
114  phi[i][p].slice(1) = dxidy_map[p]*phi_ref.slice(0) + detady_map[p]*phi_ref.slice(1)
115  + dzetady_map[p]*phi_ref.slice(2);
116 
117  phi[i][p].slice(2) = dxidz_map[p]*phi_ref.slice(0) + detadz_map[p]*phi_ref.slice(1)
118  + dzetadz_map[p]*phi_ref.slice(2);
119  }
120  break;
121  }
122 
123  default:
124  libmesh_error();
125 
126  } // switch(dim)
127 
128  return;
129  }
130 
131  template< typename OutputShape >
133  const Elem* const,
134  const std::vector<Point>&,
135  const FEGenericBase<OutputShape>& fe,
136  std::vector<std::vector<OutputShape> >& curl_phi ) const
137  {
138  switch(dim)
139  {
140  // These element transformations only make sense in 2D and 3D
141  case 0:
142  case 1:
143  {
144  libmesh_error();
145  }
146 
147  case 2:
148  {
149  const std::vector<std::vector<OutputShape> >& dphi_dxi = fe.get_dphidxi();
150  const std::vector<std::vector<OutputShape> >& dphi_deta = fe.get_dphideta();
151 
152  const std::vector<Real>& J = fe.get_fe_map().get_jacobian();
153 
154  // FIXME: I don't think this is valid for 2D elements in 3D space
155  /* In 2D: curl(phi) = J^{-1} * curl(\hat{phi}) */
156  for (unsigned int i=0; i<curl_phi.size(); i++)
157  for (unsigned int p=0; p<curl_phi[i].size(); p++)
158  {
159  curl_phi[i][p].slice(0) = curl_phi[i][p].slice(1) = 0.0;
160 
161  curl_phi[i][p].slice(2) = ( dphi_dxi[i][p].slice(1) - dphi_deta[i][p].slice(0) )/J[p];
162  }
163 
164  break;
165  }
166 
167  case 3:
168  {
169  const std::vector<std::vector<OutputShape> >& dphi_dxi = fe.get_dphidxi();
170  const std::vector<std::vector<OutputShape> >& dphi_deta = fe.get_dphideta();
171  const std::vector<std::vector<OutputShape> >& dphi_dzeta = fe.get_dphidzeta();
172 
173  const std::vector<RealGradient>& dxyz_dxi = fe.get_fe_map().get_dxyzdxi();
174  const std::vector<RealGradient>& dxyz_deta = fe.get_fe_map().get_dxyzdeta();
175  const std::vector<RealGradient>& dxyz_dzeta = fe.get_fe_map().get_dxyzdzeta();
176 
177  const std::vector<Real>& J = fe.get_fe_map().get_jacobian();
178 
179  for (unsigned int i=0; i<curl_phi.size(); i++)
180  for (unsigned int p=0; p<curl_phi[i].size(); p++)
181  {
182  Real dx_dxi = dxyz_dxi[p](0);
183  Real dx_deta = dxyz_deta[p](0);
184  Real dx_dzeta = dxyz_dzeta[p](0);
185 
186  Real dy_dxi = dxyz_dxi[p](1);
187  Real dy_deta = dxyz_deta[p](1);
188  Real dy_dzeta = dxyz_dzeta[p](1);
189 
190  Real dz_dxi = dxyz_dxi[p](2);
191  Real dz_deta = dxyz_deta[p](2);
192  Real dz_dzeta = dxyz_dzeta[p](2);
193 
194  const Real inv_jac = 1.0/J[p];
195 
196  /* In 3D: curl(phi) = J^{-1} dx/dxi * curl(\hat{phi})
197 
198  dx/dxi = [ dx/dxi dx/deta dx/dzeta
199  dy/dxi dy/deta dy/dzeta
200  dz/dxi dz/deta dz/dzeta ]
201 
202  curl(u) = [ du_z/deta - du_y/dzeta
203  du_x/dzeta - du_z/dxi
204  du_y/dxi - du_x/deta ]
205  */
206  curl_phi[i][p].slice(0) = inv_jac*( dx_dxi*( dphi_deta[i][p].slice(2) -
207  dphi_dzeta[i][p].slice(1) ) +
208  dx_deta*( dphi_dzeta[i][p].slice(0) -
209  dphi_dxi[i][p].slice(2) ) +
210  dx_dzeta*( dphi_dxi[i][p].slice(1) -
211  dphi_deta[i][p].slice(0) ) );
212 
213  curl_phi[i][p].slice(1) = inv_jac*( dy_dxi*( dphi_deta[i][p].slice(2) -
214  dphi_dzeta[i][p].slice(1) ) +
215  dy_deta*( dphi_dzeta[i][p].slice(0)-
216  dphi_dxi[i][p].slice(2) ) +
217  dy_dzeta*( dphi_dxi[i][p].slice(1) -
218  dphi_deta[i][p].slice(0) ) );
219 
220  curl_phi[i][p].slice(2) = inv_jac*( dz_dxi*( dphi_deta[i][p].slice(2) -
221  dphi_dzeta[i][p].slice(1) ) +
222  dz_deta*( dphi_dzeta[i][p].slice(0) -
223  dphi_dxi[i][p].slice(2) ) +
224  dz_dzeta*( dphi_dxi[i][p].slice(1) -
225  dphi_deta[i][p].slice(0) ) );
226  }
227 
228  break;
229  }
230 
231  default:
232  libmesh_error();
233 
234  } // switch(dim)
235 
236  return;
237  }
238 
240 
241 template<>
242 void HCurlFETransformation<Real>::map_phi( const unsigned int,
243  const Elem* const,
244  const std::vector<Point>&,
245  const FEGenericBase<Real>&,
246  std::vector<std::vector<Real> >& ) const
247 {
248  libMesh::err << "HCurl transformations only make sense for vector-valued elements."
249  << std::endl;
250  libmesh_error();
251 }
252 
253 template<>
254 void HCurlFETransformation<Real>::map_curl( const unsigned int,
255  const Elem* const,
256  const std::vector<Point>&,
257  const FEGenericBase<Real>&,
258  std::vector<std::vector<Real> >& ) const
259 {
260  libMesh::err << "HCurl transformations only make sense for vector-valued elements."
261  << std::endl;
262  libmesh_error();
263 }
264 
265 
266 } // namespace libMesh

Site Created By: libMesh Developers
Last modified: February 07 2014 16:57:05 UTC

Hosted By:
SourceForge.net Logo