libMesh
hcurl_fe_transformation.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2017 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 
18 #include "libmesh/hcurl_fe_transformation.h"
19 #include "libmesh/fe_interface.h"
20 
21 namespace libMesh
22 {
23 
24 template<typename OutputShape>
26 {
27  fe.get_fe_map().get_dxidx();
28 }
29 
30 
31 
32 template<typename OutputShape>
34 {
35  fe.get_fe_map().get_dxidx();
36 }
37 
38 
39 
40 template<typename OutputShape>
42 {
43  fe.get_fe_map().get_dxidx();
44 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
46 #endif
47 }
48 
49 
50 
51 template<typename OutputShape>
53  const Elem * const elem,
54  const std::vector<Point> & qp,
55  const FEGenericBase<OutputShape> & fe,
56  std::vector<std::vector<OutputShape>> & phi) const
57 {
58  switch (dim)
59  {
60  case 0:
61  case 1:
62  libmesh_error_msg("These element transformations only make sense in 2D and 3D.");
63 
64  case 2:
65  {
66  const std::vector<Real> & dxidx_map = fe.get_fe_map().get_dxidx();
67  const std::vector<Real> & dxidy_map = fe.get_fe_map().get_dxidy();
68 
69  const std::vector<Real> & detadx_map = fe.get_fe_map().get_detadx();
70  const std::vector<Real> & detady_map = fe.get_fe_map().get_detady();
71 
72  // FIXME: Need to update for 2D elements in 3D space
73  /* phi = (dx/dxi)^-T * \hat{phi}
74  In 2D:
75  (dx/dxi)^{-1} = [ dxi/dx dxi/dy
76  deta/dx deta/dy ]
77 
78  so: dxi/dx^{-T} * \hat{phi} = [ dxi/dx deta/dx [ \hat{phi}_xi
79  dxi/dy deta/dy ] \hat{phi}_eta ]
80 
81  or in indicial notation: phi_j = xi_{i,j}*\hat{phi}_i */
82 
83  for (std::size_t i=0; i<phi.size(); i++)
84  for (std::size_t p=0; p<phi[i].size(); p++)
85  {
86  // Need to temporarily cache reference shape functions
87  // TODO: PB: Might be worth trying to build phi_ref separately to see
88  // if we can get vectorization
89  OutputShape phi_ref;
90  FEInterface::shape<OutputShape>(2, fe.get_fe_type(), elem, i, qp[p], phi_ref);
91 
92  phi[i][p](0) = dxidx_map[p]*phi_ref.slice(0) + detadx_map[p]*phi_ref.slice(1);
93 
94  phi[i][p](1) = dxidy_map[p]*phi_ref.slice(0) + detady_map[p]*phi_ref.slice(1);
95  }
96 
97  break;
98  }
99 
100  case 3:
101  {
102  const std::vector<Real> & dxidx_map = fe.get_fe_map().get_dxidx();
103  const std::vector<Real> & dxidy_map = fe.get_fe_map().get_dxidy();
104  const std::vector<Real> & dxidz_map = fe.get_fe_map().get_dxidz();
105 
106  const std::vector<Real> & detadx_map = fe.get_fe_map().get_detadx();
107  const std::vector<Real> & detady_map = fe.get_fe_map().get_detady();
108  const std::vector<Real> & detadz_map = fe.get_fe_map().get_detadz();
109 
110  const std::vector<Real> & dzetadx_map = fe.get_fe_map().get_dzetadx();
111  const std::vector<Real> & dzetady_map = fe.get_fe_map().get_dzetady();
112  const std::vector<Real> & dzetadz_map = fe.get_fe_map().get_dzetadz();
113 
114  /* phi = (dx/dxi)^-T * \hat{phi}
115  In 3D:
116  dx/dxi^-1 = [ dxi/dx dxi/dy dxi/dz
117  deta/dx deta/dy deta/dz
118  dzeta/dx dzeta/dy dzeta/dz]
119 
120  so: dxi/dx^-T * \hat{phi} = [ dxi/dx deta/dx dzeta/dx [ \hat{phi}_xi
121  dxi/dy deta/dy dzeta/dy \hat{phi}_eta
122  dxi/dz deta/dz dzeta/dz ] \hat{phi}_zeta ]
123 
124  or in indicial notation: phi_j = xi_{i,j}*\hat{phi}_i */
125 
126  for (std::size_t i=0; i<phi.size(); i++)
127  for (std::size_t p=0; p<phi[i].size(); p++)
128  {
129  // Need to temporarily cache reference shape functions
130  // TODO: PB: Might be worth trying to build phi_ref separately to see
131  // if we can get vectorization
132  OutputShape phi_ref;
133  FEInterface::shape<OutputShape>(3, fe.get_fe_type(), elem, i, qp[p], phi_ref);
134 
135  phi[i][p].slice(0) = dxidx_map[p]*phi_ref.slice(0) + detadx_map[p]*phi_ref.slice(1)
136  + dzetadx_map[p]*phi_ref.slice(2);
137 
138  phi[i][p].slice(1) = dxidy_map[p]*phi_ref.slice(0) + detady_map[p]*phi_ref.slice(1)
139  + dzetady_map[p]*phi_ref.slice(2);
140 
141  phi[i][p].slice(2) = dxidz_map[p]*phi_ref.slice(0) + detadz_map[p]*phi_ref.slice(1)
142  + dzetadz_map[p]*phi_ref.slice(2);
143  }
144  break;
145  }
146 
147  default:
148  libmesh_error_msg("Invalid dim = " << dim);
149  } // switch(dim)
150 }
151 
152 template<typename OutputShape>
154  const Elem * const,
155  const std::vector<Point> &,
156  const FEGenericBase<OutputShape> & fe,
157  std::vector<std::vector<OutputShape>> & curl_phi) const
158 {
159  switch (dim)
160  {
161  case 0:
162  case 1:
163  libmesh_error_msg("These element transformations only make sense in 2D and 3D.");
164 
165  case 2:
166  {
167  const std::vector<std::vector<OutputShape>> & dphi_dxi = fe.get_dphidxi();
168  const std::vector<std::vector<OutputShape>> & dphi_deta = fe.get_dphideta();
169 
170  const std::vector<Real> & J = fe.get_fe_map().get_jacobian();
171 
172  // FIXME: I don't think this is valid for 2D elements in 3D space
173  /* In 2D: curl(phi) = J^{-1} * curl(\hat{phi}) */
174  for (std::size_t i=0; i<curl_phi.size(); i++)
175  for (std::size_t p=0; p<curl_phi[i].size(); p++)
176  {
177  curl_phi[i][p].slice(0) = curl_phi[i][p].slice(1) = 0.0;
178 
179  curl_phi[i][p].slice(2) = ( dphi_dxi[i][p].slice(1) - dphi_deta[i][p].slice(0) )/J[p];
180  }
181 
182  break;
183  }
184 
185  case 3:
186  {
187  const std::vector<std::vector<OutputShape>> & dphi_dxi = fe.get_dphidxi();
188  const std::vector<std::vector<OutputShape>> & dphi_deta = fe.get_dphideta();
189  const std::vector<std::vector<OutputShape>> & dphi_dzeta = fe.get_dphidzeta();
190 
191  const std::vector<RealGradient> & dxyz_dxi = fe.get_fe_map().get_dxyzdxi();
192  const std::vector<RealGradient> & dxyz_deta = fe.get_fe_map().get_dxyzdeta();
193  const std::vector<RealGradient> & dxyz_dzeta = fe.get_fe_map().get_dxyzdzeta();
194 
195  const std::vector<Real> & J = fe.get_fe_map().get_jacobian();
196 
197  for (std::size_t i=0; i<curl_phi.size(); i++)
198  for (std::size_t p=0; p<curl_phi[i].size(); p++)
199  {
200  Real dx_dxi = dxyz_dxi[p](0);
201  Real dx_deta = dxyz_deta[p](0);
202  Real dx_dzeta = dxyz_dzeta[p](0);
203 
204  Real dy_dxi = dxyz_dxi[p](1);
205  Real dy_deta = dxyz_deta[p](1);
206  Real dy_dzeta = dxyz_dzeta[p](1);
207 
208  Real dz_dxi = dxyz_dxi[p](2);
209  Real dz_deta = dxyz_deta[p](2);
210  Real dz_dzeta = dxyz_dzeta[p](2);
211 
212  const Real inv_jac = 1.0/J[p];
213 
214  /* In 3D: curl(phi) = J^{-1} dx/dxi * curl(\hat{phi})
215 
216  dx/dxi = [ dx/dxi dx/deta dx/dzeta
217  dy/dxi dy/deta dy/dzeta
218  dz/dxi dz/deta dz/dzeta ]
219 
220  curl(u) = [ du_z/deta - du_y/dzeta
221  du_x/dzeta - du_z/dxi
222  du_y/dxi - du_x/deta ]
223  */
224  curl_phi[i][p].slice(0) = inv_jac*( dx_dxi*( dphi_deta[i][p].slice(2) -
225  dphi_dzeta[i][p].slice(1) ) +
226  dx_deta*( dphi_dzeta[i][p].slice(0) -
227  dphi_dxi[i][p].slice(2) ) +
228  dx_dzeta*( dphi_dxi[i][p].slice(1) -
229  dphi_deta[i][p].slice(0) ) );
230 
231  curl_phi[i][p].slice(1) = inv_jac*( dy_dxi*( dphi_deta[i][p].slice(2) -
232  dphi_dzeta[i][p].slice(1) ) +
233  dy_deta*( dphi_dzeta[i][p].slice(0)-
234  dphi_dxi[i][p].slice(2) ) +
235  dy_dzeta*( dphi_dxi[i][p].slice(1) -
236  dphi_deta[i][p].slice(0) ) );
237 
238  curl_phi[i][p].slice(2) = inv_jac*( dz_dxi*( dphi_deta[i][p].slice(2) -
239  dphi_dzeta[i][p].slice(1) ) +
240  dz_deta*( dphi_dzeta[i][p].slice(0) -
241  dphi_dxi[i][p].slice(2) ) +
242  dz_dzeta*( dphi_dxi[i][p].slice(1) -
243  dphi_deta[i][p].slice(0) ) );
244  }
245 
246  break;
247  }
248 
249  default:
250  libmesh_error_msg("Invalid dim = " << dim);
251  } // switch(dim)
252 }
253 
255 
256 template<>
258 {
259  libmesh_error_msg("HCurl transformations only make sense for vector-valued elements.");
260 }
261 
262 template<>
264 {
265  libmesh_error_msg("HCurl transformations only make sense for vector-valued elements.");
266 }
267 
268 template<>
270 {
271  libmesh_error_msg("HCurl transformations only make sense for vector-valued elements.");
272 }
273 
274 template<>
275 void HCurlFETransformation<Real>::map_phi(const unsigned int,
276  const Elem * const,
277  const std::vector<Point> &,
278  const FEGenericBase<Real> &,
279  std::vector<std::vector<Real>> &) const
280 {
281  libmesh_error_msg("HCurl transformations only make sense for vector-valued elements.");
282 }
283 
284 template<>
286  const Elem * const,
287  const std::vector<Point> &,
288  const FEGenericBase<Real> &,
289  std::vector<std::vector<Real>> &) const
290 {
291  libmesh_error_msg("HCurl transformations only make sense for vector-valued elements.");
292 }
293 
294 
295 } // namespace libMesh
const std::vector< std::vector< Real > > & get_d2xidxyz2() const
Second derivatives of "xi" reference coordinate wrt physical coordinates.
Definition: fe_map.h:307
const std::vector< std::vector< OutputShape > > & get_dphidzeta() const
Definition: fe_base.h:280
This class handles the computation of the shape functions in the physical domain for HCurl conforming...
const std::vector< std::vector< OutputShape > > & get_dphideta() const
Definition: fe_base.h:272
unsigned int dim
const FEMap & get_fe_map() const
Definition: fe_abstract.h:452
This is the base class from which all geometric element types are derived.
Definition: elem.h:89
const std::vector< Real > & get_dxidz() const
Definition: fe_map.h:251
The libMesh namespace provides an interface to certain functionality in the library.
virtual void init_map_phi(const FEGenericBase< OutputShape > &fe) const libmesh_override
Pre-requests any necessary data from FEMap.
FEType get_fe_type() const
Definition: fe_abstract.h:421
virtual void init_map_d2phi(const FEGenericBase< OutputShape > &fe) const libmesh_override
Pre-requests any necessary data from FEMap.
const std::vector< Real > & get_dxidy() const
Definition: fe_map.h:243
const std::vector< RealGradient > & get_dxyzdxi() const
Definition: fe_map.h:161
virtual void map_curl(const unsigned int dim, const Elem *const elem, const std::vector< Point > &qp, const FEGenericBase< OutputShape > &fe, std::vector< std::vector< OutputShape >> &curl_phi) const libmesh_override
Evaluates the shape function curl in physical coordinates based on conforming finite element transfo...
const std::vector< Real > & get_jacobian() const
Definition: fe_map.h:145
virtual void map_phi(const unsigned int dim, const Elem *const elem, const std::vector< Point > &qp, const FEGenericBase< OutputShape > &fe, std::vector< std::vector< OutputShape >> &phi) const libmesh_override
Evaluates shape functions in physical coordinates for conforming elements.
const std::vector< std::vector< OutputShape > > & get_dphidxi() const
Definition: fe_base.h:264
const std::vector< RealGradient > & get_dxyzdzeta() const
Definition: fe_map.h:177
virtual void init_map_dphi(const FEGenericBase< OutputShape > &fe) const libmesh_override
Pre-requests any necessary data from FEMap.
const std::vector< Real > & get_dzetadz() const
Definition: fe_map.h:299
const std::vector< Real > & get_detadz() const
Definition: fe_map.h:275
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
const std::vector< Real > & get_dzetady() const
Definition: fe_map.h:291
const std::vector< RealGradient > & get_dxyzdeta() const
Definition: fe_map.h:169
const std::vector< Real > & get_detady() const
Definition: fe_map.h:267
const std::vector< Real > & get_dxidx() const
Definition: fe_map.h:235
const std::vector< Real > & get_detadx() const
Definition: fe_map.h:259
const std::vector< Real > & get_dzetadx() const
Definition: fe_map.h:283
This class forms the foundation from which generic finite elements may be derived.