libMesh
fe.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 
19 
20 // Local includes
21 #include "libmesh/elem.h"
22 #include "libmesh/fe.h"
23 #include "libmesh/fe_interface.h"
24 #include "libmesh/fe_macro.h"
25 #include "libmesh/libmesh_logging.h"
26 #include "libmesh/quadrature.h"
27 #include "libmesh/tensor_value.h"
28 
29 namespace libMesh
30 {
31 
32 
33 // ------------------------------------------------------------
34 // FE class members
35 template <unsigned int Dim, FEFamily T>
36 unsigned int FE<Dim,T>::n_shape_functions () const
37 {
38  return FE<Dim,T>::n_dofs (this->elem_type,
39  static_cast<Order>(this->fe_type.order + this->_p_level));
40 }
41 
42 
43 template <unsigned int Dim, FEFamily T>
45 {
46  libmesh_assert(q);
47  this->qrule = q;
48  // make sure we don't cache results from a previous quadrature rule
49  this->elem_type = INVALID_ELEM;
50  return;
51 }
52 
53 
54 template <unsigned int Dim, FEFamily T>
55 unsigned int FE<Dim,T>::n_quadrature_points () const
56 {
57  libmesh_assert(this->qrule);
58  return this->qrule->n_points();
59 }
60 
61 
62 template <unsigned int Dim, FEFamily T>
63 void FE<Dim,T>::dofs_on_side(const Elem * const elem,
64  const Order o,
65  unsigned int s,
66  std::vector<unsigned int> & di)
67 {
68  libmesh_assert(elem);
69  libmesh_assert_less (s, elem->n_sides());
70 
71  di.clear();
72  unsigned int nodenum = 0;
73  const unsigned int n_nodes = elem->n_nodes();
74  for (unsigned int n = 0; n != n_nodes; ++n)
75  {
76  const unsigned int n_dofs = n_dofs_at_node(elem->type(),
77  static_cast<Order>(o + elem->p_level()), n);
78  if (elem->is_node_on_side(n, s))
79  for (unsigned int i = 0; i != n_dofs; ++i)
80  di.push_back(nodenum++);
81  else
82  nodenum += n_dofs;
83  }
84 }
85 
86 
87 
88 template <unsigned int Dim, FEFamily T>
89 void FE<Dim,T>::dofs_on_edge(const Elem * const elem,
90  const Order o,
91  unsigned int e,
92  std::vector<unsigned int> & di)
93 {
94  libmesh_assert(elem);
95  libmesh_assert_less (e, elem->n_edges());
96 
97  di.clear();
98  unsigned int nodenum = 0;
99  const unsigned int n_nodes = elem->n_nodes();
100  for (unsigned int n = 0; n != n_nodes; ++n)
101  {
102  const unsigned int n_dofs = n_dofs_at_node(elem->type(),
103  static_cast<Order>(o + elem->p_level()), n);
104  if (elem->is_node_on_edge(n, e))
105  for (unsigned int i = 0; i != n_dofs; ++i)
106  di.push_back(nodenum++);
107  else
108  nodenum += n_dofs;
109  }
110 }
111 
112 
113 
114 template <unsigned int Dim, FEFamily T>
115 void FE<Dim,T>::reinit(const Elem * elem,
116  const std::vector<Point> * const pts,
117  const std::vector<Real> * const weights)
118 {
119  // We can be called with no element. If we're evaluating SCALAR
120  // dofs we'll still have work to do.
121  // libmesh_assert(elem);
122 
123  // We're calculating now! Time to determine what.
124  this->determine_calculations();
125 
126  // Try to avoid calling init_shape_functions
127  // even when shapes_need_reinit
128  bool cached_nodes_still_fit = false;
129 
130  // Most of the hard work happens when we have an actual element
131  if (elem)
132  {
133  // Initialize the shape functions at the user-specified
134  // points
135  if (pts != libmesh_nullptr)
136  {
137  // Set the type and p level for this element
138  this->elem_type = elem->type();
139  this->_p_level = elem->p_level();
140 
141  // Initialize the shape functions
142  this->_fe_map->template init_reference_to_physical_map<Dim>
143  (*pts, elem);
144  this->init_shape_functions (*pts, elem);
145 
146  // The shape functions do not correspond to the qrule
147  this->shapes_on_quadrature = false;
148  }
149 
150  // If there are no user specified points, we use the
151  // quadrature rule
152 
153  // update the type in accordance to the current cell
154  // and reinit if the cell type has changed or (as in
155  // the case of the hierarchics) the shape functions need
156  // reinit, since they depend on the particular element shape
157  else
158  {
159  libmesh_assert(this->qrule);
160  this->qrule->init(elem->type(), elem->p_level());
161 
162  if (this->qrule->shapes_need_reinit())
163  this->shapes_on_quadrature = false;
164 
165  if (this->elem_type != elem->type() ||
166  this->_p_level != elem->p_level() ||
167  !this->shapes_on_quadrature)
168  {
169  // Set the type and p level for this element
170  this->elem_type = elem->type();
171  this->_p_level = elem->p_level();
172  // Initialize the shape functions
173  this->_fe_map->template init_reference_to_physical_map<Dim>
174  (this->qrule->get_points(), elem);
175  this->init_shape_functions (this->qrule->get_points(), elem);
176 
177  if (this->shapes_need_reinit())
178  {
179  cached_nodes.resize(elem->n_nodes());
180  for (unsigned int n = 0; n != elem->n_nodes(); ++n)
181  {
182  cached_nodes[n] = elem->point(n);
183  }
184  }
185  }
186  else
187  {
188  // libmesh_assert_greater (elem->n_nodes(), 1);
189 
190  cached_nodes_still_fit = true;
191  if (cached_nodes.size() != elem->n_nodes())
192  cached_nodes_still_fit = false;
193  else
194  for (unsigned int n = 1; n < elem->n_nodes(); ++n)
195  {
196  if (!(elem->point(n) - elem->point(0)).relative_fuzzy_equals
197  ((cached_nodes[n] - cached_nodes[0]), 1e-13))
198  {
199  cached_nodes_still_fit = false;
200  break;
201  }
202  }
203 
204  if (this->shapes_need_reinit() && !cached_nodes_still_fit)
205  {
206  this->_fe_map->template init_reference_to_physical_map<Dim>
207  (this->qrule->get_points(), elem);
208  this->init_shape_functions (this->qrule->get_points(), elem);
209  cached_nodes.resize(elem->n_nodes());
210  for (unsigned int n = 0; n != elem->n_nodes(); ++n)
211  cached_nodes[n] = elem->point(n);
212  }
213  }
214 
215  // The shape functions correspond to the qrule
216  this->shapes_on_quadrature = true;
217  }
218  }
219  else // With no defined elem, so mapping or caching to
220  // be done, and our "quadrature rule" is one point for nonlocal
221  // (SCALAR) variables and zero points for local variables.
222  {
223  this->elem_type = INVALID_ELEM;
224  this->_p_level = 0;
225 
226  if (!pts)
227  {
228  if (T == SCALAR)
229  {
230  this->qrule->get_points() =
231  std::vector<Point>(1,Point(0));
232 
233  this->qrule->get_weights() =
234  std::vector<Real>(1,1);
235  }
236  else
237  {
238  this->qrule->get_points().clear();
239  this->qrule->get_weights().clear();
240  }
241 
242  this->init_shape_functions (this->qrule->get_points(), elem);
243  }
244  else
245  this->init_shape_functions (*pts, elem);
246  }
247 
248  // Compute the map for this element.
249  if (pts != libmesh_nullptr)
250  {
251  if (weights != libmesh_nullptr)
252  {
253  this->_fe_map->compute_map (this->dim, *weights, elem, this->calculate_d2phi);
254  }
255  else
256  {
257  std::vector<Real> dummy_weights (pts->size(), 1.);
258  this->_fe_map->compute_map (this->dim, dummy_weights, elem, this->calculate_d2phi);
259  }
260  }
261  else
262  {
263  this->_fe_map->compute_map (this->dim, this->qrule->get_weights(), elem, this->calculate_d2phi);
264  }
265 
266  // Compute the shape functions and the derivatives at all of the
267  // quadrature points.
268  if (!cached_nodes_still_fit)
269  {
270  if (pts != libmesh_nullptr)
271  this->compute_shape_functions (elem,*pts);
272  else
273  this->compute_shape_functions(elem,this->qrule->get_points());
274  }
275 }
276 
277 
278 
279 template <unsigned int Dim, FEFamily T>
280 void FE<Dim,T>::init_shape_functions(const std::vector<Point> & qp,
281  const Elem * elem)
282 {
283  // Start logging the shape function initialization
284  LOG_SCOPE("init_shape_functions()", "FE");
285 
286  // The number of quadrature points.
287  const unsigned int n_qp = cast_int<unsigned int>(qp.size());
288 
289  // Number of shape functions in the finite element approximation
290  // space.
291  const unsigned int n_approx_shape_functions =
292  this->n_shape_functions(this->get_type(),
293  this->get_order());
294 
295  // resize the vectors to hold current data
296  // Phi are the shape functions used for the FE approximation
297  // Phi_map are the shape functions used for the FE mapping
298  if (this->calculate_phi)
299  this->phi.resize (n_approx_shape_functions);
300 
301  if (this->calculate_dphi)
302  {
303  this->dphi.resize (n_approx_shape_functions);
304  this->dphidx.resize (n_approx_shape_functions);
305  this->dphidy.resize (n_approx_shape_functions);
306  this->dphidz.resize (n_approx_shape_functions);
307  }
308 
309  if (this->calculate_dphiref)
310  {
311  if (Dim > 0)
312  this->dphidxi.resize (n_approx_shape_functions);
313 
314  if (Dim > 1)
315  this->dphideta.resize (n_approx_shape_functions);
316 
317  if (Dim > 2)
318  this->dphidzeta.resize (n_approx_shape_functions);
319  }
320 
321  if (this->calculate_curl_phi && (FEInterface::field_type(T) == TYPE_VECTOR))
322  this->curl_phi.resize(n_approx_shape_functions);
323 
324  if (this->calculate_div_phi && (FEInterface::field_type(T) == TYPE_VECTOR))
325  this->div_phi.resize(n_approx_shape_functions);
326 
327 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
328  if (this->calculate_d2phi)
329  {
330  this->d2phi.resize (n_approx_shape_functions);
331  this->d2phidx2.resize (n_approx_shape_functions);
332  this->d2phidxdy.resize (n_approx_shape_functions);
333  this->d2phidxdz.resize (n_approx_shape_functions);
334  this->d2phidy2.resize (n_approx_shape_functions);
335  this->d2phidydz.resize (n_approx_shape_functions);
336  this->d2phidz2.resize (n_approx_shape_functions);
337 
338  if (Dim > 0)
339  this->d2phidxi2.resize (n_approx_shape_functions);
340 
341  if (Dim > 1)
342  {
343  this->d2phidxideta.resize (n_approx_shape_functions);
344  this->d2phideta2.resize (n_approx_shape_functions);
345  }
346  if (Dim > 2)
347  {
348  this->d2phidxidzeta.resize (n_approx_shape_functions);
349  this->d2phidetadzeta.resize (n_approx_shape_functions);
350  this->d2phidzeta2.resize (n_approx_shape_functions);
351  }
352  }
353 #endif // ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
354 
355  for (unsigned int i=0; i<n_approx_shape_functions; i++)
356  {
357  if (this->calculate_phi)
358  this->phi[i].resize (n_qp);
359  if (this->calculate_dphi)
360  {
361  this->dphi[i].resize (n_qp);
362  this->dphidx[i].resize (n_qp);
363  this->dphidy[i].resize (n_qp);
364  this->dphidz[i].resize (n_qp);
365  }
366 
367  if (this->calculate_dphiref)
368  {
369  if (Dim > 0)
370  this->dphidxi[i].resize(n_qp);
371 
372  if (Dim > 1)
373  this->dphideta[i].resize(n_qp);
374 
375  if (Dim > 2)
376  this->dphidzeta[i].resize(n_qp);
377  }
378 
379  if (this->calculate_curl_phi && (FEInterface::field_type(T) == TYPE_VECTOR))
380  this->curl_phi[i].resize(n_qp);
381 
382  if (this->calculate_div_phi && (FEInterface::field_type(T) == TYPE_VECTOR))
383  this->div_phi[i].resize(n_qp);
384 
385 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
386  if (this->calculate_d2phi)
387  {
388  this->d2phi[i].resize (n_qp);
389  this->d2phidx2[i].resize (n_qp);
390  this->d2phidxdy[i].resize (n_qp);
391  this->d2phidxdz[i].resize (n_qp);
392  this->d2phidy2[i].resize (n_qp);
393  this->d2phidydz[i].resize (n_qp);
394  this->d2phidz2[i].resize (n_qp);
395  if (Dim > 0)
396  this->d2phidxi2[i].resize (n_qp);
397  if (Dim > 1)
398  {
399  this->d2phidxideta[i].resize (n_qp);
400  this->d2phideta2[i].resize (n_qp);
401  }
402  if (Dim > 2)
403  {
404  this->d2phidxidzeta[i].resize (n_qp);
405  this->d2phidetadzeta[i].resize (n_qp);
406  this->d2phidzeta2[i].resize (n_qp);
407  }
408  }
409 #endif // ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
410  }
411 
412 
413 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
414  //------------------------------------------------------------
415  // Initialize the data fields, which should only be used for infinite
416  // elements, to some sensible values, so that using a FE with the
417  // variational formulation of an InfFE, correct element matrices are
418  // returned
419 
420  {
421  this->weight.resize (n_qp);
422  this->dweight.resize (n_qp);
423  this->dphase.resize (n_qp);
424 
425  for (unsigned int p=0; p<n_qp; p++)
426  {
427  this->weight[p] = 1.;
428  this->dweight[p].zero();
429  this->dphase[p].zero();
430  }
431 
432  }
433 #endif // ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
434 
435  switch (Dim)
436  {
437 
438  //------------------------------------------------------------
439  // 0D
440  case 0:
441  {
442  break;
443  }
444 
445  //------------------------------------------------------------
446  // 1D
447  case 1:
448  {
449  // Compute the value of the approximation shape function i at quadrature point p
450  if (this->calculate_dphiref)
451  for (unsigned int i=0; i<n_approx_shape_functions; i++)
452  for (unsigned int p=0; p<n_qp; p++)
453  this->dphidxi[i][p] = FE<Dim,T>::shape_deriv (elem, this->fe_type.order, i, 0, qp[p]);
454 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
455  if (this->calculate_d2phi)
456  for (unsigned int i=0; i<n_approx_shape_functions; i++)
457  for (unsigned int p=0; p<n_qp; p++)
458  this->d2phidxi2[i][p] = FE<Dim,T>::shape_second_deriv (elem, this->fe_type.order, i, 0, qp[p]);
459 #endif // ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
460 
461  break;
462  }
463 
464 
465 
466  //------------------------------------------------------------
467  // 2D
468  case 2:
469  {
470  // Compute the value of the approximation shape function i at quadrature point p
471  if (this->calculate_dphiref)
472  for (unsigned int i=0; i<n_approx_shape_functions; i++)
473  for (unsigned int p=0; p<n_qp; p++)
474  {
475  this->dphidxi[i][p] = FE<Dim,T>::shape_deriv (elem, this->fe_type.order, i, 0, qp[p]);
476  this->dphideta[i][p] = FE<Dim,T>::shape_deriv (elem, this->fe_type.order, i, 1, qp[p]);
477  }
478 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
479  if (this->calculate_d2phi)
480  for (unsigned int i=0; i<n_approx_shape_functions; i++)
481  for (unsigned int p=0; p<n_qp; p++)
482  {
483  this->d2phidxi2[i][p] = FE<Dim,T>::shape_second_deriv (elem, this->fe_type.order, i, 0, qp[p]);
484  this->d2phidxideta[i][p] = FE<Dim,T>::shape_second_deriv (elem, this->fe_type.order, i, 1, qp[p]);
485  this->d2phideta2[i][p] = FE<Dim,T>::shape_second_deriv (elem, this->fe_type.order, i, 2, qp[p]);
486  }
487 #endif // ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
488 
489 
490  break;
491  }
492 
493 
494 
495  //------------------------------------------------------------
496  // 3D
497  case 3:
498  {
499  // Compute the value of the approximation shape function i at quadrature point p
500  if (this->calculate_dphiref)
501  for (unsigned int i=0; i<n_approx_shape_functions; i++)
502  for (unsigned int p=0; p<n_qp; p++)
503  {
504  this->dphidxi[i][p] = FE<Dim,T>::shape_deriv (elem, this->fe_type.order, i, 0, qp[p]);
505  this->dphideta[i][p] = FE<Dim,T>::shape_deriv (elem, this->fe_type.order, i, 1, qp[p]);
506  this->dphidzeta[i][p] = FE<Dim,T>::shape_deriv (elem, this->fe_type.order, i, 2, qp[p]);
507  }
508 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
509  if (this->calculate_d2phi)
510  for (unsigned int i=0; i<n_approx_shape_functions; i++)
511  for (unsigned int p=0; p<n_qp; p++)
512  {
513  this->d2phidxi2[i][p] = FE<Dim,T>::shape_second_deriv (elem, this->fe_type.order, i, 0, qp[p]);
514  this->d2phidxideta[i][p] = FE<Dim,T>::shape_second_deriv (elem, this->fe_type.order, i, 1, qp[p]);
515  this->d2phideta2[i][p] = FE<Dim,T>::shape_second_deriv (elem, this->fe_type.order, i, 2, qp[p]);
516  this->d2phidxidzeta[i][p] = FE<Dim,T>::shape_second_deriv (elem, this->fe_type.order, i, 3, qp[p]);
517  this->d2phidetadzeta[i][p] = FE<Dim,T>::shape_second_deriv (elem, this->fe_type.order, i, 4, qp[p]);
518  this->d2phidzeta2[i][p] = FE<Dim,T>::shape_second_deriv (elem, this->fe_type.order, i, 5, qp[p]);
519  }
520 #endif // ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
521 
522  break;
523  }
524 
525 
526  default:
527  libmesh_error_msg("Invalid dimension Dim = " << Dim);
528  }
529 }
530 
531 
532 
533 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
534 
535 template <unsigned int Dim, FEFamily T>
536 void FE<Dim,T>::init_base_shape_functions(const std::vector<Point> & qp,
537  const Elem * e)
538 {
539  this->elem_type = e->type();
540  this->_fe_map->template init_reference_to_physical_map<Dim>(qp, e);
541  init_shape_functions(qp, e);
542 }
543 
544 #endif // LIBMESH_ENABLE_INFINITE_ELEMENTS
545 
546 
547 
548 //--------------------------------------------------------------
549 // Explicit instantiations using macro from fe_macro.h
550 
551 INSTANTIATE_FE(0);
552 
553 INSTANTIATE_FE(1);
554 
555 INSTANTIATE_FE(2);
556 
557 INSTANTIATE_FE(3);
558 
560 
561 } // namespace libMesh
virtual bool is_node_on_edge(const unsigned int n, const unsigned int e) const =0
virtual bool is_node_on_side(const unsigned int n, const unsigned int s) const =0
static unsigned int n_dofs(const ElemType t, const Order o)
virtual ElemType type() const =0
unsigned int p_level() const
Definition: elem.h:2422
INSTANTIATE_SUBDIVISION_FE
Definition: fe.C:559
unsigned int dim
virtual unsigned int n_edges() const =0
This is the base class from which all geometric element types are derived.
Definition: elem.h:89
static FEFieldType field_type(const FEType &fe_type)
const class libmesh_nullptr_t libmesh_nullptr
static OutputShape shape_deriv(const ElemType t, const Order o, const unsigned int i, const unsigned int j, const Point &p)
static void dofs_on_side(const Elem *const elem, const Order o, unsigned int s, std::vector< unsigned int > &di)
Fills the vector di with the local degree of freedom indices associated with side s of element elem...
Definition: fe.C:63
The libMesh namespace provides an interface to certain functionality in the library.
libmesh_assert(j)
virtual unsigned int n_nodes() const =0
A specific instantiation of the FEBase class.
Definition: fe.h:89
const dof_id_type n_nodes
Definition: tecplot_io.C:67
virtual unsigned int n_quadrature_points() const libmesh_override
Definition: fe.C:55
INSTANTIATE_FE(0)
dof_id_type weight(const MeshBase &mesh, const processor_id_type pid)
Definition: mesh_tools.C:245
virtual void attach_quadrature_rule(QBase *q) libmesh_override
Provides the class with the quadrature rule, which provides the locations (on a reference element) wh...
Definition: fe.C:44
virtual unsigned int n_sides() const =0
static void dofs_on_edge(const Elem *const elem, const Order o, unsigned int e, std::vector< unsigned int > &di)
Fills the vector di with the local degree of freedom indices associated with edge e of element elem...
Definition: fe.C:89
const Point & point(const unsigned int i) const
Definition: elem.h:1809
virtual void reinit(const Elem *elem, const std::vector< Point > *const pts=libmesh_nullptr, const std::vector< Real > *const weights=libmesh_nullptr) libmesh_override
This is at the core of this class.
Definition: fe.C:115
virtual unsigned int n_shape_functions() const libmesh_override
Definition: fe.C:36
virtual void init_shape_functions(const std::vector< Point > &qp, const Elem *e)
Update the various member data fields phi, dphidxi, dphideta, dphidzeta, etc.
Definition: fe.C:280
virtual void init_base_shape_functions(const std::vector< Point > &qp, const Elem *e) libmesh_override
Initialize the data fields for the base of an an infinite element.
Definition: fe.C:536
Order
defines an enum for polynomial orders.
Definition: enum_order.h:32
A Point defines a location in LIBMESH_DIM dimensional Real space.
Definition: point.h:38
The QBase class provides the basic functionality from which various quadrature rules can be derived...
Definition: quadrature.h:53
static OutputShape shape_second_deriv(const ElemType t, const Order o, const unsigned int i, const unsigned int j, const Point &p)