fe.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 
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"
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  libmesh_assert(elem);
120 
121  // Try to avoid calling init_shape_functions
122  // even when shapes_need_reinit
123  bool cached_nodes_still_fit = false;
124 
125  // Initialize the shape functions at the user-specified
126  // points
127  if (pts != NULL)
128  {
129  // Set the type and p level for this element
130  this->elem_type = elem->type();
131  this->_p_level = elem->p_level();
132 
133  // Initialize the shape functions
134  this->_fe_map->template init_reference_to_physical_map<Dim>(*pts, elem);
135  this->init_shape_functions (*pts, elem);
136 
137  // The shape functions do not correspond to the qrule
138  this->shapes_on_quadrature = false;
139  }
140 
141  // If there are no user specified points, we use the
142  // quadrature rule
143 
144  // update the type in accordance to the current cell
145  // and reinit if the cell type has changed or (as in
146  // the case of the hierarchics) the shape functions need
147  // reinit, since they depend on the particular element shape
148  else
149  {
150  libmesh_assert(this->qrule);
151  this->qrule->init(elem->type(), elem->p_level());
152 
153  if(this->qrule->shapes_need_reinit())
154  this->shapes_on_quadrature = false;
155 
156  if (this->elem_type != elem->type() ||
157  this->_p_level != elem->p_level() ||
158  !this->shapes_on_quadrature)
159  {
160  // Set the type and p level for this element
161  this->elem_type = elem->type();
162  this->_p_level = elem->p_level();
163  // Initialize the shape functions
164  this->_fe_map->template init_reference_to_physical_map<Dim>(this->qrule->get_points(), elem);
165  this->init_shape_functions (this->qrule->get_points(), elem);
166 
167  if (this->shapes_need_reinit())
168  {
169  cached_nodes.resize(elem->n_nodes());
170  for (unsigned int n = 0; n != elem->n_nodes(); ++n)
171  {
172  cached_nodes[n] = elem->point(n);
173  }
174  }
175  }
176  else
177  {
178  // libmesh_assert_greater (elem->n_nodes(), 1);
179 
180  cached_nodes_still_fit = true;
181  if (cached_nodes.size() != elem->n_nodes())
182  cached_nodes_still_fit = false;
183  else
184  for (unsigned int n = 1; n < elem->n_nodes(); ++n)
185  {
186  if (!(elem->point(n) - elem->point(0)).relative_fuzzy_equals(
187  (cached_nodes[n] - cached_nodes[0]), 1e-13))
188  {
189  cached_nodes_still_fit = false;
190  break;
191  }
192  }
193 
194  if (this->shapes_need_reinit() && !cached_nodes_still_fit)
195  {
196  this->_fe_map->template init_reference_to_physical_map<Dim>(this->qrule->get_points(), elem);
197  this->init_shape_functions (this->qrule->get_points(), elem);
198  cached_nodes.resize(elem->n_nodes());
199  for (unsigned int n = 0; n != elem->n_nodes(); ++n)
200  cached_nodes[n] = elem->point(n);
201  }
202  }
203 
204  // The shape functions correspond to the qrule
205  this->shapes_on_quadrature = true;
206  }
207 
208  // Compute the map for this element. In the future we can specify
209  // different types of maps
210  if (pts != NULL)
211  {
212  if (weights != NULL)
213  {
214  this->_fe_map->compute_map (this->dim,*weights, elem);
215  }
216  else
217  {
218  std::vector<Real> dummy_weights (pts->size(), 1.);
219  this->_fe_map->compute_map (this->dim,dummy_weights, elem);
220  }
221  }
222  else
223  {
224  this->_fe_map->compute_map (this->dim,this->qrule->get_weights(), elem);
225  }
226 
227  // Compute the shape functions and the derivatives at all of the
228  // quadrature points.
229  if (!cached_nodes_still_fit)
230  {
231  if (pts != NULL)
232  this->compute_shape_functions (elem,*pts);
233  else
234  this->compute_shape_functions(elem,this->qrule->get_points());
235  }
236 }
237 
238 
239 
240 template <unsigned int Dim, FEFamily T>
241 void FE<Dim,T>::init_shape_functions(const std::vector<Point>& qp,
242  const Elem* elem)
243 {
244  libmesh_assert(elem);
245  this->calculations_started = true;
246 
247  // If the user forgot to request anything, we'll be safe and
248  // calculate everything:
249 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
250  if (!this->calculate_phi && !this->calculate_dphi && !this->calculate_d2phi
251  && !this->calculate_curl_phi && !this->calculate_div_phi)
252  {
253  this->calculate_phi = this->calculate_dphi = this->calculate_d2phi = this->calculate_dphiref = true;
255  {
256  this->calculate_curl_phi = true;
257  this->calculate_div_phi = true;
258  }
259  }
260 #else
261  if (!this->calculate_phi && !this->calculate_dphi && !this->calculate_curl_phi && !this->calculate_div_phi)
262  {
263  this->calculate_phi = this->calculate_dphi = this->calculate_dphiref = true;
265  {
266  this->calculate_curl_phi = true;
267  this->calculate_div_phi = true;
268  }
269  }
270 #endif // LIBMESH_ENABLE_SECOND_DERIVATIVES
271 
272  // Start logging the shape function initialization
273  START_LOG("init_shape_functions()", "FE");
274 
275 
276  // The number of quadrature points.
277  const unsigned int n_qp = libmesh_cast_int<unsigned int>(qp.size());
278 
279  // Number of shape functions in the finite element approximation
280  // space.
281  const unsigned int n_approx_shape_functions =
282  this->n_shape_functions(this->get_type(),
283  this->get_order());
284 
285  // resize the vectors to hold current data
286  // Phi are the shape functions used for the FE approximation
287  // Phi_map are the shape functions used for the FE mapping
288  if (this->calculate_phi)
289  this->phi.resize (n_approx_shape_functions);
290 
291  if (this->calculate_dphi)
292  {
293  this->dphi.resize (n_approx_shape_functions);
294  this->dphidx.resize (n_approx_shape_functions);
295  this->dphidy.resize (n_approx_shape_functions);
296  this->dphidz.resize (n_approx_shape_functions);
297  }
298 
299  if(this->calculate_dphiref)
300  {
301  if (Dim > 0)
302  this->dphidxi.resize (n_approx_shape_functions);
303 
304  if (Dim > 1)
305  this->dphideta.resize (n_approx_shape_functions);
306 
307  if (Dim > 2)
308  this->dphidzeta.resize (n_approx_shape_functions);
309  }
310 
311  if( this->calculate_curl_phi && (FEInterface::field_type(T) == TYPE_VECTOR) )
312  this->curl_phi.resize(n_approx_shape_functions);
313 
314  if( this->calculate_div_phi && (FEInterface::field_type(T) == TYPE_VECTOR) )
315  this->div_phi.resize(n_approx_shape_functions);
316 
317 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
318  if (this->calculate_d2phi)
319  {
320  this->d2phi.resize (n_approx_shape_functions);
321  this->d2phidx2.resize (n_approx_shape_functions);
322  this->d2phidxdy.resize (n_approx_shape_functions);
323  this->d2phidxdz.resize (n_approx_shape_functions);
324  this->d2phidy2.resize (n_approx_shape_functions);
325  this->d2phidydz.resize (n_approx_shape_functions);
326  this->d2phidz2.resize (n_approx_shape_functions);
327 
328  if (Dim > 0)
329  this->d2phidxi2.resize (n_approx_shape_functions);
330 
331  if (Dim > 1)
332  {
333  this->d2phidxideta.resize (n_approx_shape_functions);
334  this->d2phideta2.resize (n_approx_shape_functions);
335  }
336  if (Dim > 2)
337  {
338  this->d2phidxidzeta.resize (n_approx_shape_functions);
339  this->d2phidetadzeta.resize (n_approx_shape_functions);
340  this->d2phidzeta2.resize (n_approx_shape_functions);
341  }
342  }
343 #endif // ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
344 
345  for (unsigned int i=0; i<n_approx_shape_functions; i++)
346  {
347  if (this->calculate_phi)
348  this->phi[i].resize (n_qp);
349  if (this->calculate_dphi)
350  {
351  this->dphi[i].resize (n_qp);
352  this->dphidx[i].resize (n_qp);
353  this->dphidy[i].resize (n_qp);
354  this->dphidz[i].resize (n_qp);
355  }
356 
357  if(this->calculate_dphiref)
358  {
359  if (Dim > 0)
360  this->dphidxi[i].resize(n_qp);
361 
362  if (Dim > 1)
363  this->dphideta[i].resize(n_qp);
364 
365  if (Dim > 2)
366  this->dphidzeta[i].resize(n_qp);
367  }
368 
369  if(this->calculate_curl_phi && (FEInterface::field_type(T) == TYPE_VECTOR) )
370  this->curl_phi[i].resize(n_qp);
371 
372  if(this->calculate_div_phi && (FEInterface::field_type(T) == TYPE_VECTOR) )
373  this->div_phi[i].resize(n_qp);
374 
375 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
376  if (this->calculate_d2phi)
377  {
378  this->d2phi[i].resize (n_qp);
379  this->d2phidx2[i].resize (n_qp);
380  this->d2phidxdy[i].resize (n_qp);
381  this->d2phidxdz[i].resize (n_qp);
382  this->d2phidy2[i].resize (n_qp);
383  this->d2phidydz[i].resize (n_qp);
384  this->d2phidz2[i].resize (n_qp);
385  if (Dim > 0)
386  this->d2phidxi2[i].resize (n_qp);
387  if (Dim > 1)
388  {
389  this->d2phidxideta[i].resize (n_qp);
390  this->d2phideta2[i].resize (n_qp);
391  }
392  if (Dim > 2)
393  {
394  this->d2phidxidzeta[i].resize (n_qp);
395  this->d2phidetadzeta[i].resize (n_qp);
396  this->d2phidzeta2[i].resize (n_qp);
397  }
398  }
399 #endif // ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
400  }
401 
402 
403 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
404  //------------------------------------------------------------
405  // Initialize the data fields, which should only be used for infinite
406  // elements, to some sensible values, so that using a FE with the
407  // variational formulation of an InfFE, correct element matrices are
408  // returned
409 
410  {
411  this->weight.resize (n_qp);
412  this->dweight.resize (n_qp);
413  this->dphase.resize (n_qp);
414 
415  for (unsigned int p=0; p<n_qp; p++)
416  {
417  this->weight[p] = 1.;
418  this->dweight[p].zero();
419  this->dphase[p].zero();
420  }
421 
422  }
423 #endif // ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
424 
425  switch (Dim)
426  {
427 
428  //------------------------------------------------------------
429  // 0D
430  case 0:
431  {
432  break;
433  }
434 
435  //------------------------------------------------------------
436  // 1D
437  case 1:
438  {
439  // Compute the value of the approximation shape function i at quadrature point p
440  if (this->calculate_dphiref)
441  for (unsigned int i=0; i<n_approx_shape_functions; i++)
442  for (unsigned int p=0; p<n_qp; p++)
443  this->dphidxi[i][p] = FE<Dim,T>::shape_deriv (elem, this->fe_type.order, i, 0, qp[p]);
444 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
445  if (this->calculate_d2phi)
446  for (unsigned int i=0; i<n_approx_shape_functions; i++)
447  for (unsigned int p=0; p<n_qp; p++)
448  this->d2phidxi2[i][p] = FE<Dim,T>::shape_second_deriv (elem, this->fe_type.order, i, 0, qp[p]);
449 #endif // ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
450 
451  break;
452  }
453 
454 
455 
456  //------------------------------------------------------------
457  // 2D
458  case 2:
459  {
460  // Compute the value of the approximation shape function i at quadrature point p
461  if (this->calculate_dphiref)
462  for (unsigned int i=0; i<n_approx_shape_functions; i++)
463  for (unsigned int p=0; p<n_qp; p++)
464  {
465  this->dphidxi[i][p] = FE<Dim,T>::shape_deriv (elem, this->fe_type.order, i, 0, qp[p]);
466  this->dphideta[i][p] = FE<Dim,T>::shape_deriv (elem, this->fe_type.order, i, 1, qp[p]);
467  }
468 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
469  if (this->calculate_d2phi)
470  for (unsigned int i=0; i<n_approx_shape_functions; i++)
471  for (unsigned int p=0; p<n_qp; p++)
472  {
473  this->d2phidxi2[i][p] = FE<Dim,T>::shape_second_deriv (elem, this->fe_type.order, i, 0, qp[p]);
474  this->d2phidxideta[i][p] = FE<Dim,T>::shape_second_deriv (elem, this->fe_type.order, i, 1, qp[p]);
475  this->d2phideta2[i][p] = FE<Dim,T>::shape_second_deriv (elem, this->fe_type.order, i, 2, qp[p]);
476  }
477 #endif // ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
478 
479 
480  break;
481  }
482 
483 
484 
485  //------------------------------------------------------------
486  // 3D
487  case 3:
488  {
489  // Compute the value of the approximation shape function i at quadrature point p
490  if (this->calculate_dphiref)
491  for (unsigned int i=0; i<n_approx_shape_functions; i++)
492  for (unsigned int p=0; p<n_qp; p++)
493  {
494  this->dphidxi[i][p] = FE<Dim,T>::shape_deriv (elem, this->fe_type.order, i, 0, qp[p]);
495  this->dphideta[i][p] = FE<Dim,T>::shape_deriv (elem, this->fe_type.order, i, 1, qp[p]);
496  this->dphidzeta[i][p] = FE<Dim,T>::shape_deriv (elem, this->fe_type.order, i, 2, qp[p]);
497  }
498 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
499  if (this->calculate_d2phi)
500  for (unsigned int i=0; i<n_approx_shape_functions; i++)
501  for (unsigned int p=0; p<n_qp; p++)
502  {
503  this->d2phidxi2[i][p] = FE<Dim,T>::shape_second_deriv (elem, this->fe_type.order, i, 0, qp[p]);
504  this->d2phidxideta[i][p] = FE<Dim,T>::shape_second_deriv (elem, this->fe_type.order, i, 1, qp[p]);
505  this->d2phideta2[i][p] = FE<Dim,T>::shape_second_deriv (elem, this->fe_type.order, i, 2, qp[p]);
506  this->d2phidxidzeta[i][p] = FE<Dim,T>::shape_second_deriv (elem, this->fe_type.order, i, 3, qp[p]);
507  this->d2phidetadzeta[i][p] = FE<Dim,T>::shape_second_deriv (elem, this->fe_type.order, i, 4, qp[p]);
508  this->d2phidzeta2[i][p] = FE<Dim,T>::shape_second_deriv (elem, this->fe_type.order, i, 5, qp[p]);
509  }
510 #endif // ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
511 
512  break;
513  }
514 
515 
516  default:
517  libmesh_error();
518  }
519 
520  // Stop logging the shape function initialization
521  STOP_LOG("init_shape_functions()", "FE");
522 }
523 
524 
525 
526 
527 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
528 
529 template <unsigned int Dim, FEFamily T>
530 void FE<Dim,T>::init_base_shape_functions(const std::vector<Point>& qp,
531  const Elem* e)
532 {
533  // I don't understand infinite elements well enough to risk
534  // calculating too little. :-( RHS
535  this->calculate_phi = this->calculate_dphi = this->calculate_d2phi = true;
536 
537  this->elem_type = e->type();
538  this->_fe_map->template init_reference_to_physical_map<Dim>(qp, e);
539  init_shape_functions(qp, e);
540 }
541 
542 #endif // LIBMESH_ENABLE_INFINITE_ELEMENTS
543 
544 
545 
546 //--------------------------------------------------------------
547 // Explicit instantiations using macro from fe_macro.h
548 
549 INSTANTIATE_FE(0);
550 
551 INSTANTIATE_FE(1);
552 
553 INSTANTIATE_FE(2);
554 
555 INSTANTIATE_FE(3);
556 
557 
558 } // namespace libMesh

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

Hosted By:
SourceForge.net Logo