libMesh
Loading...
Searching...
No Matches
hcurl_fe_transformation.C
Go to the documentation of this file.
1// The libMesh Finite Element Library.
2// Copyright (C) 2002-2026 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
22namespace libMesh
23{
24
25template<typename OutputShape>
30
31
32
33template<typename OutputShape>
38
39
40
41template<typename OutputShape>
42void 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
52template<typename OutputShape>
53void 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,
58 const bool /*add_p_level*/) const
59{
60 switch (dim)
61 {
62 case 0:
63 case 1:
64 libmesh_error_msg("These element transformations only make sense in 2D and 3D.");
65
66 case 2:
67 {
68 const std::vector<Real> & dxidx_map = fe.get_fe_map().get_dxidx();
69 const std::vector<Real> & dxidy_map = fe.get_fe_map().get_dxidy();
70 const std::vector<Real> & dxidz_map = fe.get_fe_map().get_dxidz();
71
72 const std::vector<Real> & detadx_map = fe.get_fe_map().get_detadx();
73 const std::vector<Real> & detady_map = fe.get_fe_map().get_detady();
74 const std::vector<Real> & detadz_map = fe.get_fe_map().get_detadz();
75
76 // phi = (dx/dxi)^-T * \hat{phi}
77 // In 2D (embedded in 3D):
78 // (dx/dxi)^{-1} = [ dxi/dx dxi/dy dxi/dz
79 // deta/dx deta/dy deta/dz ]
80 //
81 // so: dxi/dx^{-T} * \hat{phi} = [ dxi/dx deta/dx [ \hat{phi}_xi
82 // dxi/dy deta/dy \hat{phi}_eta ]
83 // dxi/dz deta/dz ]
84 //
85 // or in indicial notation: phi_j = xi_{i,j}*\hat{phi}_i
86 for (auto i : index_range(phi))
87 for (auto p : index_range(phi[i]))
88 {
89 // Need to temporarily cache reference shape functions
90 // TODO: PB: Might be worth trying to build phi_ref separately to see
91 // if we can get vectorization
92 // We are computing mapping basis functions, so we explicitly ignore
93 // any non-zero p_level() the Elem might have.
94 OutputShape phi_ref;
95 FEInterface::shape(fe.get_fe_type(), /*extra_order=*/0, elem, i, qp[p], phi_ref);
96
97 phi[i][p](0) = dxidx_map[p]*phi_ref.slice(0) + detadx_map[p]*phi_ref.slice(1);
98
99 phi[i][p](1) = dxidy_map[p]*phi_ref.slice(0) + detady_map[p]*phi_ref.slice(1);
100
101 phi[i][p](2) = dxidz_map[p]*phi_ref.slice(0) + detadz_map[p]*phi_ref.slice(1);
102 }
103
104 break;
105 }
106
107 case 3:
108 {
109 const std::vector<Real> & dxidx_map = fe.get_fe_map().get_dxidx();
110 const std::vector<Real> & dxidy_map = fe.get_fe_map().get_dxidy();
111 const std::vector<Real> & dxidz_map = fe.get_fe_map().get_dxidz();
112
113 const std::vector<Real> & detadx_map = fe.get_fe_map().get_detadx();
114 const std::vector<Real> & detady_map = fe.get_fe_map().get_detady();
115 const std::vector<Real> & detadz_map = fe.get_fe_map().get_detadz();
116
117 const std::vector<Real> & dzetadx_map = fe.get_fe_map().get_dzetadx();
118 const std::vector<Real> & dzetady_map = fe.get_fe_map().get_dzetady();
119 const std::vector<Real> & dzetadz_map = fe.get_fe_map().get_dzetadz();
120
121 // phi = (dx/dxi)^-T * \hat{phi}
122 // In 3D:
123 // dx/dxi^-1 = [ dxi/dx dxi/dy dxi/dz
124 // deta/dx deta/dy deta/dz
125 // dzeta/dx dzeta/dy dzeta/dz]
126 //
127 // so: dxi/dx^-T * \hat{phi} = [ dxi/dx deta/dx dzeta/dx [ \hat{phi}_xi
128 // dxi/dy deta/dy dzeta/dy \hat{phi}_eta
129 // dxi/dz deta/dz dzeta/dz ] \hat{phi}_zeta ]
130 //
131 // or in indicial notation: phi_j = xi_{i,j}*\hat{phi}_i
132 for (auto i : index_range(phi))
133 for (auto p : index_range(phi[i]))
134 {
135 // Need to temporarily cache reference shape functions
136 // TODO: PB: Might be worth trying to build phi_ref separately to see
137 // if we can get vectorization
138 // We are computing mapping basis functions, so we explicitly ignore
139 // any non-zero p_level() the Elem might have.
140 OutputShape phi_ref;
141 FEInterface::shape(fe.get_fe_type(), /*extra_order=*/0, elem, i, qp[p], phi_ref);
142
143 phi[i][p].slice(0) = dxidx_map[p]*phi_ref.slice(0) + detadx_map[p]*phi_ref.slice(1)
144 + dzetadx_map[p]*phi_ref.slice(2);
145
146 phi[i][p].slice(1) = dxidy_map[p]*phi_ref.slice(0) + detady_map[p]*phi_ref.slice(1)
147 + dzetady_map[p]*phi_ref.slice(2);
148
149 phi[i][p].slice(2) = dxidz_map[p]*phi_ref.slice(0) + detadz_map[p]*phi_ref.slice(1)
150 + dzetadz_map[p]*phi_ref.slice(2);
151 }
152 break;
153 }
154
155 default:
156 libmesh_error_msg("Invalid dim = " << dim);
157 } // switch(dim)
158}
159
160template<typename OutputShape>
161void HCurlFETransformation<OutputShape>::map_curl(const unsigned int dim,
162 const Elem * const,
163 const std::vector<Point> &,
164 const FEGenericBase<OutputShape> & fe,
165 std::vector<std::vector<OutputShape>> & curl_phi) const
166{
167 switch (dim)
168 {
169 case 0:
170 case 1:
171 libmesh_error_msg("These element transformations only make sense in 2D and 3D.");
172
173 case 2:
174 {
175 const std::vector<Real> & dxidx_map = fe.get_fe_map().get_dxidx();
176 const std::vector<Real> & dxidy_map = fe.get_fe_map().get_dxidy();
177 const std::vector<Real> & dxidz_map = fe.get_fe_map().get_dxidz();
178
179 const std::vector<Real> & detadx_map = fe.get_fe_map().get_detadx();
180 const std::vector<Real> & detady_map = fe.get_fe_map().get_detady();
181 const std::vector<Real> & detadz_map = fe.get_fe_map().get_detadz();
182
183 const std::vector<std::vector<OutputShape>> & dphi_dxi = fe.get_dphidxi();
184 const std::vector<std::vector<OutputShape>> & dphi_deta = fe.get_dphideta();
185
186 // In 2D (embedded in 3D): curl(phi)_j = C(j) * curl(\hat{phi}) where C(j)
187 // is the jth column cofactor of (dx/dxi)^{-1} = [ dxi/dx dxi/dy dxi/dz
188 // deta/dx deta/dy deta/dz ]
189 for (auto i : index_range(curl_phi))
190 for (auto p : index_range(curl_phi[i]))
191 {
192 const Real curl_ref = dphi_dxi[i][p].slice(1) - dphi_deta[i][p].slice(0);
193 curl_phi[i][p].slice(0) = curl_ref * (dxidy_map[p]*detadz_map[p] - dxidz_map[p]*detady_map[p]);
194 curl_phi[i][p].slice(1) = curl_ref * -(dxidx_map[p]*detadz_map[p] - dxidz_map[p]*detadx_map[p]);
195 curl_phi[i][p].slice(2) = curl_ref * (dxidx_map[p]*detady_map[p] - dxidy_map[p]*detadx_map[p]);
196 }
197
198 break;
199 }
200
201 case 3:
202 {
203 const std::vector<std::vector<OutputShape>> & dphi_dxi = fe.get_dphidxi();
204 const std::vector<std::vector<OutputShape>> & dphi_deta = fe.get_dphideta();
205 const std::vector<std::vector<OutputShape>> & dphi_dzeta = fe.get_dphidzeta();
206
207 const std::vector<RealGradient> & dxyz_dxi = fe.get_fe_map().get_dxyzdxi();
208 const std::vector<RealGradient> & dxyz_deta = fe.get_fe_map().get_dxyzdeta();
209 const std::vector<RealGradient> & dxyz_dzeta = fe.get_fe_map().get_dxyzdzeta();
210
211 const std::vector<Real> & J = fe.get_fe_map().get_jacobian();
212
213 for (auto i : index_range(curl_phi))
214 for (auto p : index_range(curl_phi[i]))
215 {
216 Real dx_dxi = dxyz_dxi[p](0);
217 Real dx_deta = dxyz_deta[p](0);
218 Real dx_dzeta = dxyz_dzeta[p](0);
219
220 Real dy_dxi = dxyz_dxi[p](1);
221 Real dy_deta = dxyz_deta[p](1);
222 Real dy_dzeta = dxyz_dzeta[p](1);
223
224 Real dz_dxi = dxyz_dxi[p](2);
225 Real dz_deta = dxyz_deta[p](2);
226 Real dz_dzeta = dxyz_dzeta[p](2);
227
228 const Real inv_jac = 1.0/J[p];
229
230 /* In 3D: curl(phi) = J^{-1} dx/dxi * curl(\hat{phi})
231
232 dx/dxi = [ dx/dxi dx/deta dx/dzeta
233 dy/dxi dy/deta dy/dzeta
234 dz/dxi dz/deta dz/dzeta ]
235
236 curl(u) = [ du_z/deta - du_y/dzeta
237 du_x/dzeta - du_z/dxi
238 du_y/dxi - du_x/deta ]
239 */
240 curl_phi[i][p].slice(0) = inv_jac*( dx_dxi*( dphi_deta[i][p].slice(2) -
241 dphi_dzeta[i][p].slice(1) ) +
242 dx_deta*( dphi_dzeta[i][p].slice(0) -
243 dphi_dxi[i][p].slice(2) ) +
244 dx_dzeta*( dphi_dxi[i][p].slice(1) -
245 dphi_deta[i][p].slice(0) ) );
246
247 curl_phi[i][p].slice(1) = inv_jac*( dy_dxi*( dphi_deta[i][p].slice(2) -
248 dphi_dzeta[i][p].slice(1) ) +
249 dy_deta*( dphi_dzeta[i][p].slice(0)-
250 dphi_dxi[i][p].slice(2) ) +
251 dy_dzeta*( dphi_dxi[i][p].slice(1) -
252 dphi_deta[i][p].slice(0) ) );
253
254 curl_phi[i][p].slice(2) = inv_jac*( dz_dxi*( dphi_deta[i][p].slice(2) -
255 dphi_dzeta[i][p].slice(1) ) +
256 dz_deta*( dphi_dzeta[i][p].slice(0) -
257 dphi_dxi[i][p].slice(2) ) +
258 dz_dzeta*( dphi_dxi[i][p].slice(1) -
259 dphi_deta[i][p].slice(0) ) );
260 }
261
262 break;
263 }
264
265 default:
266 libmesh_error_msg("Invalid dim = " << dim);
267 } // switch(dim)
268}
269
270template class LIBMESH_EXPORT HCurlFETransformation<RealGradient>;
271
272template<>
273void HCurlFETransformation<Real>::init_map_phi(const FEGenericBase<Real> & ) const
274{
275 libmesh_error_msg("HCurl transformations only make sense for vector-valued elements.");
276}
277
278template<>
279void HCurlFETransformation<Real>::init_map_dphi(const FEGenericBase<Real> & ) const
280{
281 libmesh_error_msg("HCurl transformations only make sense for vector-valued elements.");
282}
283
284template<>
285void HCurlFETransformation<Real>::init_map_d2phi(const FEGenericBase<Real> & ) const
286{
287 libmesh_error_msg("HCurl transformations only make sense for vector-valued elements.");
288}
289
290template<>
291void HCurlFETransformation<Real>::map_phi(const unsigned int,
292 const Elem * const,
293 const std::vector<Point> &,
294 const FEGenericBase<Real> &,
295 std::vector<std::vector<Real>> &,
296 bool) const
297{
298 libmesh_error_msg("HCurl transformations only make sense for vector-valued elements.");
299}
300
301template<>
302void HCurlFETransformation<Real>::map_curl(const unsigned int,
303 const Elem * const,
304 const std::vector<Point> &,
305 const FEGenericBase<Real> &,
306 std::vector<std::vector<Real>> &) const
307{
308 libmesh_error_msg("HCurl transformations only make sense for vector-valued elements.");
309}
310
311
312} // namespace libMesh
dim
unsigned int dim
Definition adaptivity_ex3.C:114
libMesh::Elem
This is the base class from which all geometric element types are derived.
Definition elem.h:96
libMesh::FEAbstract::get_fe_type
FEType get_fe_type() const
Definition fe_abstract.h:520
libMesh::FEAbstract::get_fe_map
const FEMap & get_fe_map() const
Definition fe_abstract.h:554
libMesh::FEGenericBase
This class forms the foundation from which generic finite elements may be derived.
Definition fe_base.h:86
libMesh::FEGenericBase::get_dphidxi
const std::vector< std::vector< OutputShape > > & get_dphidxi() const
Definition fe_base.h:293
libMesh::FEGenericBase::get_dphidzeta
const std::vector< std::vector< OutputShape > > & get_dphidzeta() const
Definition fe_base.h:309
libMesh::FEGenericBase::get_dphideta
const std::vector< std::vector< OutputShape > > & get_dphideta() const
Definition fe_base.h:301
libMesh::FEInterface::shape
static Real shape(const unsigned int dim, const FEType &fe_t, const ElemType t, const unsigned int i, const Point &p)
Definition fe_interface.C:761
libMesh::FEMap::get_dxidz
const std::vector< Real > & get_dxidz() const
Definition fe_map.h:325
libMesh::FEMap::get_dxyzdeta
const std::vector< RealGradient > & get_dxyzdeta() const
Definition fe_map.h:247
libMesh::FEMap::get_detadx
const std::vector< Real > & get_detadx() const
Definition fe_map.h:333
libMesh::FEMap::get_jacobian
const std::vector< Real > & get_jacobian() const
Definition fe_map.h:223
libMesh::FEMap::get_dzetady
const std::vector< Real > & get_dzetady() const
Definition fe_map.h:365
libMesh::FEMap::get_detady
const std::vector< Real > & get_detady() const
Definition fe_map.h:341
libMesh::FEMap::get_dxidx
const std::vector< Real > & get_dxidx() const
Definition fe_map.h:309
libMesh::FEMap::get_dzetadx
const std::vector< Real > & get_dzetadx() const
Definition fe_map.h:357
libMesh::FEMap::get_dxyzdzeta
const std::vector< RealGradient > & get_dxyzdzeta() const
Definition fe_map.h:255
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:381
libMesh::FEMap::get_dxidy
const std::vector< Real > & get_dxidy() const
Definition fe_map.h:317
libMesh::FEMap::get_dxyzdxi
const std::vector< RealGradient > & get_dxyzdxi() const
Definition fe_map.h:239
libMesh::FEMap::get_dzetadz
const std::vector< Real > & get_dzetadz() const
Definition fe_map.h:373
libMesh::FEMap::get_detadz
const std::vector< Real > & get_detadz() const
Definition fe_map.h:349
libMesh::HCurlFETransformation
This class handles the computation of the shape functions in the physical domain for HCurl conforming...
Definition hcurl_fe_transformation.h:37
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::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, bool add_p_level=true) const override
Evaluates shape functions in physical coordinates for conforming elements.
Definition hcurl_fe_transformation.C:53
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:161
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::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
libMesh
The libMesh namespace provides an interface to certain functionality in the library.
libMesh::index_range
auto index_range(const T &sizable)
Helper function that returns an IntRange<std::size_t> representing all the indices of the passed-in v...
Definition int_range.h:153
libMesh::Real
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
Definition libmesh_common.h:144
Generated on Fri Aug 21 2026 13:32:29 for libMesh by doxygen 1.9.8