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