libMesh
fe_xyz.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/libmesh_logging.h"
22 #include "libmesh/fe.h"
23 #include "libmesh/elem.h"
24 #include "libmesh/fe_interface.h"
25 #include "libmesh/string_to_enum.h"
26 
27 namespace libMesh
28 {
29 
30 // ------------------------------------------------------------
31 // XYZ-specific implementations
32 
33 // Anonymous namespace for local helper functions
34 namespace {
35 
36 void xyz_nodal_soln(const Elem * elem,
37  const Order order,
38  const std::vector<Number> & elem_soln,
39  std::vector<Number> & nodal_soln,
40  unsigned Dim)
41 {
42  const unsigned int n_nodes = elem->n_nodes();
43 
44  const ElemType elem_type = elem->type();
45 
46  nodal_soln.resize(n_nodes);
47 
48  const Order totalorder = static_cast<Order>(order + elem->p_level());
49 
50  switch (totalorder)
51  {
52  // Constant shape functions
53  case CONSTANT:
54  {
55  libmesh_assert_equal_to (elem_soln.size(), 1);
56 
57  const Number val = elem_soln[0];
58 
59  for (unsigned int n=0; n<n_nodes; n++)
60  nodal_soln[n] = val;
61 
62  return;
63  }
64 
65 
66  // For other orders do interpolation at the nodes
67  // explicitly.
68  default:
69  {
70  // FEType object to be passed to various FEInterface functions below.
71  FEType fe_type(totalorder, XYZ);
72 
73  const unsigned int n_sf =
74  // FE<Dim,T>::n_shape_functions(elem_type, totalorder);
75  FEInterface::n_shape_functions(Dim, fe_type, elem_type);
76 
77  for (unsigned int n=0; n<n_nodes; n++)
78  {
79  libmesh_assert_equal_to (elem_soln.size(), n_sf);
80 
81  // Zero before summation
82  nodal_soln[n] = 0;
83 
84  // u_i = Sum (alpha_i phi_i)
85  for (unsigned int i=0; i<n_sf; i++)
86  nodal_soln[n] += elem_soln[i] *
87  // FE<Dim,T>::shape(elem, order, i, elem->point(n));
88  FEInterface::shape(Dim, fe_type, elem, i, elem->point(n));
89  }
90 
91  return;
92  } // default
93  } // switch
94 } // xyz_nodal_soln()
95 
96 
97 
98 
99 
100 unsigned int xyz_n_dofs(const ElemType t, const Order o)
101 {
102  switch (o)
103  {
104 
105  // constant shape functions
106  // no matter what shape there is only one DOF.
107  case CONSTANT:
108  return (t != INVALID_ELEM) ? 1 : 0;
109 
110 
111  // Discontinuous linear shape functions
112  // expressed in the XYZ monomials.
113  case FIRST:
114  {
115  switch (t)
116  {
117  case NODEELEM:
118  return 1;
119 
120  case EDGE2:
121  case EDGE3:
122  case EDGE4:
123  return 2;
124 
125  case TRI3:
126  case TRISHELL3:
127  case TRI6:
128  case QUAD4:
129  case QUADSHELL4:
130  case QUAD8:
131  case QUADSHELL8:
132  case QUAD9:
133  return 3;
134 
135  case TET4:
136  case TET10:
137  case HEX8:
138  case HEX20:
139  case HEX27:
140  case PRISM6:
141  case PRISM15:
142  case PRISM18:
143  case PYRAMID5:
144  case PYRAMID13:
145  case PYRAMID14:
146  return 4;
147 
148  case INVALID_ELEM:
149  return 0;
150 
151  default:
152  libmesh_error_msg("ERROR: Bad ElemType = " << Utility::enum_to_string(t) << " for " << Utility::enum_to_string(o) << " order approximation!");
153  }
154  }
155 
156 
157  // Discontinuous quadratic shape functions
158  // expressed in the XYZ monomials.
159  case SECOND:
160  {
161  switch (t)
162  {
163  case NODEELEM:
164  return 1;
165 
166  case EDGE2:
167  case EDGE3:
168  case EDGE4:
169  return 3;
170 
171  case TRI3:
172  case TRISHELL3:
173  case TRI6:
174  case QUAD4:
175  case QUADSHELL4:
176  case QUAD8:
177  case QUADSHELL8:
178  case QUAD9:
179  return 6;
180 
181  case TET4:
182  case TET10:
183  case HEX8:
184  case HEX20:
185  case HEX27:
186  case PRISM6:
187  case PRISM15:
188  case PRISM18:
189  case PYRAMID5:
190  case PYRAMID13:
191  case PYRAMID14:
192  return 10;
193 
194  case INVALID_ELEM:
195  return 0;
196 
197  default:
198  libmesh_error_msg("ERROR: Bad ElemType = " << Utility::enum_to_string(t) << " for " << Utility::enum_to_string(o) << " order approximation!");
199  }
200  }
201 
202 
203  // Discontinuous cubic shape functions
204  // expressed in the XYZ monomials.
205  case THIRD:
206  {
207  switch (t)
208  {
209  case NODEELEM:
210  return 1;
211 
212  case EDGE2:
213  case EDGE3:
214  case EDGE4:
215  return 4;
216 
217  case TRI3:
218  case TRISHELL3:
219  case TRI6:
220  case QUAD4:
221  case QUADSHELL4:
222  case QUAD8:
223  case QUADSHELL8:
224  case QUAD9:
225  return 10;
226 
227  case TET4:
228  case TET10:
229  case HEX8:
230  case HEX20:
231  case HEX27:
232  case PRISM6:
233  case PRISM15:
234  case PRISM18:
235  case PYRAMID5:
236  case PYRAMID13:
237  case PYRAMID14:
238  return 20;
239 
240  case INVALID_ELEM:
241  return 0;
242 
243  default:
244  libmesh_error_msg("ERROR: Bad ElemType = " << Utility::enum_to_string(t) << " for " << Utility::enum_to_string(o) << " order approximation!");
245  }
246  }
247 
248 
249  // Discontinuous quartic shape functions
250  // expressed in the XYZ monomials.
251  case FOURTH:
252  {
253  switch (t)
254  {
255  case NODEELEM:
256  return 1;
257 
258  case EDGE2:
259  case EDGE3:
260  return 5;
261 
262  case TRI3:
263  case TRISHELL3:
264  case TRI6:
265  case QUAD4:
266  case QUADSHELL4:
267  case QUAD8:
268  case QUADSHELL8:
269  case QUAD9:
270  return 15;
271 
272  case TET4:
273  case TET10:
274  case HEX8:
275  case HEX20:
276  case HEX27:
277  case PRISM6:
278  case PRISM15:
279  case PRISM18:
280  case PYRAMID5:
281  case PYRAMID13:
282  case PYRAMID14:
283  return 35;
284 
285  case INVALID_ELEM:
286  return 0;
287 
288  default:
289  libmesh_error_msg("ERROR: Bad ElemType = " << Utility::enum_to_string(t) << " for " << Utility::enum_to_string(o) << " order approximation!");
290  }
291  }
292 
293 
294  default:
295  {
296  const unsigned int order = static_cast<unsigned int>(o);
297  switch (t)
298  {
299  case NODEELEM:
300  return 1;
301 
302  case EDGE2:
303  case EDGE3:
304  return (order+1);
305 
306  case TRI3:
307  case TRISHELL3:
308  case TRI6:
309  case QUAD4:
310  case QUADSHELL4:
311  case QUAD8:
312  case QUADSHELL8:
313  case QUAD9:
314  return (order+1)*(order+2)/2;
315 
316  case TET4:
317  case TET10:
318  case HEX8:
319  case HEX20:
320  case HEX27:
321  case PRISM6:
322  case PRISM15:
323  case PRISM18:
324  case PYRAMID5:
325  case PYRAMID13:
326  case PYRAMID14:
327  return (order+1)*(order+2)*(order+3)/6;
328 
329  case INVALID_ELEM:
330  return 0;
331 
332  default:
333  libmesh_error_msg("ERROR: Bad ElemType = " << Utility::enum_to_string(t) << " for " << Utility::enum_to_string(o) << " order approximation!");
334  }
335  }
336  }
337 
338  libmesh_error_msg("We'll never get here!");
339  return 0;
340 }
341 
342 
343 
344 
345 unsigned int xyz_n_dofs_per_elem(const ElemType t,
346  const Order o)
347 {
348  switch (o)
349  {
350  // constant shape functions always have 1 DOF per element
351  case CONSTANT:
352  return (t != INVALID_ELEM) ? 1 : 0;
353 
354 
355  // Discontinuous linear shape functions
356  // expressed in the XYZ monomials.
357  case FIRST:
358  {
359  switch (t)
360  {
361  case NODEELEM:
362  return 1;
363 
364  // 1D linears have 2 DOFs per element
365  case EDGE2:
366  case EDGE3:
367  case EDGE4:
368  return 2;
369 
370  // 2D linears have 3 DOFs per element
371  case TRI3:
372  case TRISHELL3:
373  case TRI6:
374  case QUAD4:
375  case QUADSHELL4:
376  case QUAD8:
377  case QUADSHELL8:
378  case QUAD9:
379  return 3;
380 
381  // 3D linears have 4 DOFs per element
382  case TET4:
383  case TET10:
384  case HEX8:
385  case HEX20:
386  case HEX27:
387  case PRISM6:
388  case PRISM15:
389  case PRISM18:
390  case PYRAMID5:
391  case PYRAMID13:
392  case PYRAMID14:
393  return 4;
394 
395  case INVALID_ELEM:
396  return 0;
397 
398  default:
399  libmesh_error_msg("ERROR: Bad ElemType = " << Utility::enum_to_string(t) << " for " << Utility::enum_to_string(o) << " order approximation!");
400  }
401  }
402 
403 
404  // Discontinuous quadratic shape functions
405  // expressed in the XYZ monomials.
406  case SECOND:
407  {
408  switch (t)
409  {
410  case NODEELEM:
411  return 1;
412 
413  // 1D quadratics have 3 DOFs per element
414  case EDGE2:
415  case EDGE3:
416  case EDGE4:
417  return 3;
418 
419  // 2D quadratics have 6 DOFs per element
420  case TRI3:
421  case TRISHELL3:
422  case TRI6:
423  case QUAD4:
424  case QUADSHELL4:
425  case QUAD8:
426  case QUADSHELL8:
427  case QUAD9:
428  return 6;
429 
430  // 3D quadratics have 10 DOFs per element
431  case TET4:
432  case TET10:
433  case HEX8:
434  case HEX20:
435  case HEX27:
436  case PRISM6:
437  case PRISM15:
438  case PRISM18:
439  case PYRAMID5:
440  case PYRAMID13:
441  case PYRAMID14:
442  return 10;
443 
444  case INVALID_ELEM:
445  return 0;
446 
447  default:
448  libmesh_error_msg("ERROR: Bad ElemType = " << Utility::enum_to_string(t) << " for " << Utility::enum_to_string(o) << " order approximation!");
449  }
450  }
451 
452 
453  // Discontinuous cubic shape functions
454  // expressed in the XYZ monomials.
455  case THIRD:
456  {
457  switch (t)
458  {
459  case NODEELEM:
460  return 1;
461 
462  case EDGE2:
463  case EDGE3:
464  case EDGE4:
465  return 4;
466 
467  case TRI3:
468  case TRISHELL3:
469  case TRI6:
470  case QUAD4:
471  case QUADSHELL4:
472  case QUAD8:
473  case QUADSHELL8:
474  case QUAD9:
475  return 10;
476 
477  case TET4:
478  case TET10:
479  case HEX8:
480  case HEX20:
481  case HEX27:
482  case PRISM6:
483  case PRISM15:
484  case PRISM18:
485  case PYRAMID5:
486  case PYRAMID13:
487  case PYRAMID14:
488  return 20;
489 
490  case INVALID_ELEM:
491  return 0;
492 
493  default:
494  libmesh_error_msg("ERROR: Bad ElemType = " << Utility::enum_to_string(t) << " for " << Utility::enum_to_string(o) << " order approximation!");
495  }
496  }
497 
498 
499  // Discontinuous quartic shape functions
500  // expressed in the XYZ monomials.
501  case FOURTH:
502  {
503  switch (t)
504  {
505  case NODEELEM:
506  return 1;
507 
508  case EDGE2:
509  case EDGE3:
510  case EDGE4:
511  return 5;
512 
513  case TRI3:
514  case TRISHELL3:
515  case TRI6:
516  case QUAD4:
517  case QUADSHELL4:
518  case QUAD8:
519  case QUADSHELL8:
520  case QUAD9:
521  return 15;
522 
523  case TET4:
524  case TET10:
525  case HEX8:
526  case HEX20:
527  case HEX27:
528  case PRISM6:
529  case PRISM15:
530  case PRISM18:
531  case PYRAMID5:
532  case PYRAMID13:
533  case PYRAMID14:
534  return 35;
535 
536  case INVALID_ELEM:
537  return 0;
538 
539  default:
540  libmesh_error_msg("ERROR: Bad ElemType = " << Utility::enum_to_string(t) << " for " << Utility::enum_to_string(o) << " order approximation!");
541  }
542  }
543 
544  default:
545  {
546  const unsigned int order = static_cast<unsigned int>(o);
547  switch (t)
548  {
549  case NODEELEM:
550  return 1;
551 
552  case EDGE2:
553  case EDGE3:
554  return (order+1);
555 
556  case TRI3:
557  case TRISHELL3:
558  case TRI6:
559  case QUAD4:
560  case QUADSHELL4:
561  case QUAD8:
562  case QUADSHELL8:
563  case QUAD9:
564  return (order+1)*(order+2)/2;
565 
566  case TET4:
567  case TET10:
568  case HEX8:
569  case HEX20:
570  case HEX27:
571  case PRISM6:
572  case PRISM15:
573  case PRISM18:
574  case PYRAMID5:
575  case PYRAMID13:
576  case PYRAMID14:
577  return (order+1)*(order+2)*(order+3)/6;
578 
579  case INVALID_ELEM:
580  return 0;
581 
582  default:
583  libmesh_error_msg("ERROR: Bad ElemType = " << Utility::enum_to_string(t) << " for " << Utility::enum_to_string(o) << " order approximation!");
584  }
585  }
586  return 0;
587  }
588 }
589 
590 
591 } // anonymous namespace
592 
593 
594 
595 
596 
597 
598 
599 template <unsigned int Dim>
600 void FEXYZ<Dim>::init_shape_functions(const std::vector<Point> & qp,
601  const Elem * libmesh_dbg_var(elem))
602 {
603  libmesh_assert(elem);
604  this->calculations_started = true;
605 
606  // If the user forgot to request anything, we'll be safe and
607  // calculate everything:
608 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
609  if (!this->calculate_phi && !this->calculate_dphi && !this->calculate_d2phi)
610  this->calculate_phi = this->calculate_dphi = this->calculate_d2phi = true;
611 #else
612  if (!this->calculate_phi && !this->calculate_dphi)
613  this->calculate_phi = this->calculate_dphi = true;
614 #endif // LIBMESH_ENABLE_SECOND_DERIVATIVES
615 
616  // Start logging the shape function initialization
617  LOG_SCOPE("init_shape_functions()", "FE");
618 
619  // The number of quadrature points.
620  const std::size_t n_qp = qp.size();
621 
622  // Number of shape functions in the finite element approximation
623  // space.
624  const unsigned int n_approx_shape_functions =
625  this->n_shape_functions(this->get_type(),
626  this->get_order());
627 
628  // resize the vectors to hold current data
629  // Phi are the shape functions used for the FE approximation
630  {
631  // (note: GCC 3.4.0 requires the use of this-> here)
632  if (this->calculate_phi)
633  this->phi.resize (n_approx_shape_functions);
634  if (this->calculate_dphi)
635  {
636  this->dphi.resize (n_approx_shape_functions);
637  this->dphidx.resize (n_approx_shape_functions);
638  this->dphidy.resize (n_approx_shape_functions);
639  this->dphidz.resize (n_approx_shape_functions);
640  }
641 
642 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
643  if (this->calculate_d2phi)
644  {
645  this->d2phi.resize (n_approx_shape_functions);
646  this->d2phidx2.resize (n_approx_shape_functions);
647  this->d2phidxdy.resize (n_approx_shape_functions);
648  this->d2phidxdz.resize (n_approx_shape_functions);
649  this->d2phidy2.resize (n_approx_shape_functions);
650  this->d2phidydz.resize (n_approx_shape_functions);
651  this->d2phidz2.resize (n_approx_shape_functions);
652  this->d2phidxi2.resize (n_approx_shape_functions);
653  }
654 #endif // ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
655 
656  for (unsigned int i=0; i<n_approx_shape_functions; i++)
657  {
658  if (this->calculate_phi)
659  this->phi[i].resize (n_qp);
660  if (this->calculate_dphi)
661  {
662  this->dphi[i].resize (n_qp);
663  this->dphidx[i].resize (n_qp);
664  this->dphidy[i].resize (n_qp);
665  this->dphidz[i].resize (n_qp);
666  }
667 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
668  if (this->calculate_d2phi)
669  {
670  this->d2phi[i].resize (n_qp);
671  this->d2phidx2[i].resize (n_qp);
672  this->d2phidxdy[i].resize (n_qp);
673  this->d2phidxdz[i].resize (n_qp);
674  this->d2phidy2[i].resize (n_qp);
675  this->d2phidydz[i].resize (n_qp);
676  this->d2phidz2[i].resize (n_qp);
677  }
678 #endif // ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
679  }
680  }
681 
682 
683 
684 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
685  //------------------------------------------------------------
686  // Initialize the data fields, which should only be used for infinite
687  // elements, to some sensible values, so that using a FE with the
688  // variational formulation of an InfFE, correct element matrices are
689  // returned
690 
691  {
692  this->weight.resize (n_qp);
693  this->dweight.resize (n_qp);
694  this->dphase.resize (n_qp);
695 
696  for (unsigned int p=0; p<n_qp; p++)
697  {
698  this->weight[p] = 1.;
699  this->dweight[p].zero();
700  this->dphase[p].zero();
701  }
702 
703  }
704 #endif // ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
705 }
706 
707 
708 
709 
710 template <unsigned int Dim>
712  const std::vector<Point> &)
713 {
714  libmesh_assert(elem);
715 
716  //-------------------------------------------------------------------------
717  // Compute the shape function values (and derivatives)
718  // at the Quadrature points. Note that the actual values
719  // have already been computed via init_shape_functions
720 
721  // Start logging the shape function computation
722  LOG_SCOPE("compute_shape_functions()", "FE");
723 
724  const std::vector<Point> & xyz_qp = this->get_xyz();
725 
726  // Compute the value of the derivative shape function i at quadrature point p
727  switch (this->dim)
728  {
729 
730  case 1:
731  {
732  if (this->calculate_phi)
733  for (std::size_t i=0; i<this->phi.size(); i++)
734  for (std::size_t p=0; p<this->phi[i].size(); p++)
735  this->phi[i][p] = FE<Dim,XYZ>::shape (elem, this->fe_type.order, i, xyz_qp[p]);
736  if (this->calculate_dphi)
737  for (std::size_t i=0; i<this->dphi.size(); i++)
738  for (std::size_t p=0; p<this->dphi[i].size(); p++)
739  {
740  this->dphi[i][p](0) =
741  this->dphidx[i][p] = FE<Dim,XYZ>::shape_deriv (elem, this->fe_type.order, i, 0, xyz_qp[p]);
742 
743  this->dphi[i][p](1) = this->dphidy[i][p] = 0.;
744  this->dphi[i][p](2) = this->dphidz[i][p] = 0.;
745  }
746 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
747  if (this->calculate_d2phi)
748  for (std::size_t i=0; i<this->d2phi.size(); i++)
749  for (std::size_t p=0; p<this->d2phi[i].size(); p++)
750  {
751  this->d2phi[i][p](0,0) =
752  this->d2phidx2[i][p] = FE<Dim,XYZ>::shape_second_deriv (elem, this->fe_type.order, i, 0, xyz_qp[p]);
753 
754 #if LIBMESH_DIM>1
755  this->d2phi[i][p](0,1) = this->d2phidxdy[i][p] =
756  this->d2phi[i][p](1,0) = 0.;
757  this->d2phi[i][p](1,1) = this->d2phidy2[i][p] = 0.;
758 #if LIBMESH_DIM>2
759  this->d2phi[i][p](0,2) = this->d2phidxdz[i][p] =
760  this->d2phi[i][p](2,0) = 0.;
761  this->d2phi[i][p](1,2) = this->d2phidydz[i][p] =
762  this->d2phi[i][p](2,1) = 0.;
763  this->d2phi[i][p](2,2) = this->d2phidz2[i][p] = 0.;
764 #endif
765 #endif
766  }
767 #endif // ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
768 
769  // All done
770  break;
771  }
772 
773  case 2:
774  {
775  if (this->calculate_phi)
776  for (std::size_t i=0; i<this->phi.size(); i++)
777  for (std::size_t p=0; p<this->phi[i].size(); p++)
778  this->phi[i][p] = FE<Dim,XYZ>::shape (elem, this->fe_type.order, i, xyz_qp[p]);
779  if (this->calculate_dphi)
780  for (std::size_t i=0; i<this->dphi.size(); i++)
781  for (std::size_t p=0; p<this->dphi[i].size(); p++)
782  {
783  this->dphi[i][p](0) =
784  this->dphidx[i][p] = FE<Dim,XYZ>::shape_deriv (elem, this->fe_type.order, i, 0, xyz_qp[p]);
785 
786  this->dphi[i][p](1) =
787  this->dphidy[i][p] = FE<Dim,XYZ>::shape_deriv (elem, this->fe_type.order, i, 1, xyz_qp[p]);
788 
789 #if LIBMESH_DIM == 3
790  this->dphi[i][p](2) = // can only assign to the Z component if LIBMESH_DIM==3
791 #endif
792  this->dphidz[i][p] = 0.;
793  }
794 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
795  if (this->calculate_d2phi)
796  for (std::size_t i=0; i<this->d2phi.size(); i++)
797  for (std::size_t p=0; p<this->d2phi[i].size(); p++)
798  {
799  this->d2phi[i][p](0,0) =
800  this->d2phidx2[i][p] = FE<Dim,XYZ>::shape_second_deriv (elem, this->fe_type.order, i, 0, xyz_qp[p]);
801 
802  this->d2phi[i][p](0,1) = this->d2phidxdy[i][p] =
803  this->d2phi[i][p](1,0) = FE<Dim,XYZ>::shape_second_deriv (elem, this->fe_type.order, i, 1, xyz_qp[p]);
804  this->d2phi[i][p](1,1) =
805  this->d2phidy2[i][p] = FE<Dim,XYZ>::shape_second_deriv (elem, this->fe_type.order, i, 2, xyz_qp[p]);
806 #if LIBMESH_DIM>2
807  this->d2phi[i][p](0,2) = this->d2phidxdz[i][p] =
808  this->d2phi[i][p](2,0) = 0.;
809  this->d2phi[i][p](1,2) = this->d2phidydz[i][p] =
810  this->d2phi[i][p](2,1) = 0.;
811  this->d2phi[i][p](2,2) = this->d2phidz2[i][p] = 0.;
812 #endif
813  }
814 #endif // ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
815 
816  // All done
817  break;
818  }
819 
820  case 3:
821  {
822  if (this->calculate_phi)
823  for (std::size_t i=0; i<this->phi.size(); i++)
824  for (std::size_t p=0; p<this->phi[i].size(); p++)
825  this->phi[i][p] = FE<Dim,XYZ>::shape (elem, this->fe_type.order, i, xyz_qp[p]);
826 
827  if (this->calculate_dphi)
828  for (std::size_t i=0; i<this->dphi.size(); i++)
829  for (std::size_t p=0; p<this->dphi[i].size(); p++)
830  {
831  this->dphi[i][p](0) =
832  this->dphidx[i][p] = FE<Dim,XYZ>::shape_deriv (elem, this->fe_type.order, i, 0, xyz_qp[p]);
833 
834  this->dphi[i][p](1) =
835  this->dphidy[i][p] = FE<Dim,XYZ>::shape_deriv (elem, this->fe_type.order, i, 1, xyz_qp[p]);
836 
837  this->dphi[i][p](2) =
838  this->dphidz[i][p] = FE<Dim,XYZ>::shape_deriv (elem, this->fe_type.order, i, 2, xyz_qp[p]);
839  }
840 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
841  if (this->calculate_d2phi)
842  for (std::size_t i=0; i<this->d2phi.size(); i++)
843  for (std::size_t p=0; p<this->d2phi[i].size(); p++)
844  {
845  this->d2phi[i][p](0,0) =
846  this->d2phidx2[i][p] = FE<Dim,XYZ>::shape_second_deriv (elem, this->fe_type.order, i, 0, xyz_qp[p]);
847 
848  this->d2phi[i][p](0,1) = this->d2phidxdy[i][p] =
849  this->d2phi[i][p](1,0) = FE<Dim,XYZ>::shape_second_deriv (elem, this->fe_type.order, i, 1, xyz_qp[p]);
850  this->d2phi[i][p](1,1) =
851  this->d2phidy2[i][p] = FE<Dim,XYZ>::shape_second_deriv (elem, this->fe_type.order, i, 2, xyz_qp[p]);
852  this->d2phi[i][p](0,2) = this->d2phidxdz[i][p] =
853  this->d2phi[i][p](2,0) = FE<Dim,XYZ>::shape_second_deriv (elem, this->fe_type.order, i, 3, xyz_qp[p]);
854  this->d2phi[i][p](1,2) = this->d2phidydz[i][p] =
855  this->d2phi[i][p](2,1) = FE<Dim,XYZ>::shape_second_deriv (elem, this->fe_type.order, i, 4, xyz_qp[p]);
856  this->d2phi[i][p](2,2) = this->d2phidz2[i][p] = FE<Dim,XYZ>::shape_second_deriv (elem, this->fe_type.order, i, 5, xyz_qp[p]);
857  }
858 #endif // ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
859 
860  // All done
861  break;
862  }
863 
864  default:
865  libmesh_error_msg("ERROR: Invalid dimension " << this->dim);
866  }
867 }
868 
869 
870 
871 
872 // Do full-specialization for every dimension, instead
873 // of explicit instantiation at the end of this file.
874 // This could be macro-ified so that it fits on one line...
875 template <>
876 void FE<0,XYZ>::nodal_soln(const Elem * elem,
877  const Order order,
878  const std::vector<Number> & elem_soln,
879  std::vector<Number> & nodal_soln)
880 { xyz_nodal_soln(elem, order, elem_soln, nodal_soln, /*Dim=*/0); }
881 
882 template <>
883 void FE<1,XYZ>::nodal_soln(const Elem * elem,
884  const Order order,
885  const std::vector<Number> & elem_soln,
886  std::vector<Number> & nodal_soln)
887 { xyz_nodal_soln(elem, order, elem_soln, nodal_soln, /*Dim=*/1); }
888 
889 template <>
890 void FE<2,XYZ>::nodal_soln(const Elem * elem,
891  const Order order,
892  const std::vector<Number> & elem_soln,
893  std::vector<Number> & nodal_soln)
894 { xyz_nodal_soln(elem, order, elem_soln, nodal_soln, /*Dim=*/2); }
895 
896 template <>
897 void FE<3,XYZ>::nodal_soln(const Elem * elem,
898  const Order order,
899  const std::vector<Number> & elem_soln,
900  std::vector<Number> & nodal_soln)
901 { xyz_nodal_soln(elem, order, elem_soln, nodal_soln, /*Dim=*/3); }
902 
903 
904 
905 // Full specialization of n_dofs() function for every dimension
906 template <> unsigned int FE<0,XYZ>::n_dofs(const ElemType t, const Order o) { return xyz_n_dofs(t, o); }
907 template <> unsigned int FE<1,XYZ>::n_dofs(const ElemType t, const Order o) { return xyz_n_dofs(t, o); }
908 template <> unsigned int FE<2,XYZ>::n_dofs(const ElemType t, const Order o) { return xyz_n_dofs(t, o); }
909 template <> unsigned int FE<3,XYZ>::n_dofs(const ElemType t, const Order o) { return xyz_n_dofs(t, o); }
910 
911 // Full specialization of n_dofs_at_node() function for every dimension.
912 // XYZ FEMs have no dofs at nodes, only element dofs.
913 template <> unsigned int FE<0,XYZ>::n_dofs_at_node(const ElemType, const Order, const unsigned int) { return 0; }
914 template <> unsigned int FE<1,XYZ>::n_dofs_at_node(const ElemType, const Order, const unsigned int) { return 0; }
915 template <> unsigned int FE<2,XYZ>::n_dofs_at_node(const ElemType, const Order, const unsigned int) { return 0; }
916 template <> unsigned int FE<3,XYZ>::n_dofs_at_node(const ElemType, const Order, const unsigned int) { return 0; }
917 
918 // Full specialization of n_dofs_per_elem() function for every dimension.
919 template <> unsigned int FE<0,XYZ>::n_dofs_per_elem(const ElemType t, const Order o) { return xyz_n_dofs_per_elem(t, o); }
920 template <> unsigned int FE<1,XYZ>::n_dofs_per_elem(const ElemType t, const Order o) { return xyz_n_dofs_per_elem(t, o); }
921 template <> unsigned int FE<2,XYZ>::n_dofs_per_elem(const ElemType t, const Order o) { return xyz_n_dofs_per_elem(t, o); }
922 template <> unsigned int FE<3,XYZ>::n_dofs_per_elem(const ElemType t, const Order o) { return xyz_n_dofs_per_elem(t, o); }
923 
924 // Full specialization of get_continuity() function for every dimension.
925 template <> FEContinuity FE<0,XYZ>::get_continuity() const { return DISCONTINUOUS; }
926 template <> FEContinuity FE<1,XYZ>::get_continuity() const { return DISCONTINUOUS; }
927 template <> FEContinuity FE<2,XYZ>::get_continuity() const { return DISCONTINUOUS; }
928 template <> FEContinuity FE<3,XYZ>::get_continuity() const { return DISCONTINUOUS; }
929 
930 // Full specialization of is_hierarchic() function for every dimension.
931 // The XYZ shape functions are hierarchic!
932 template <> bool FE<0,XYZ>::is_hierarchic() const { return true; }
933 template <> bool FE<1,XYZ>::is_hierarchic() const { return true; }
934 template <> bool FE<2,XYZ>::is_hierarchic() const { return true; }
935 template <> bool FE<3,XYZ>::is_hierarchic() const { return true; }
936 
937 #ifdef LIBMESH_ENABLE_AMR
938 
939 // Full specialization of compute_constraints() function for 2D and
940 // 3D only. There are no constraints for discontinuous elements, so
941 // we do nothing.
942 template <> void FE<2,XYZ>::compute_constraints (DofConstraints &, DofMap &, const unsigned int, const Elem *) {}
943 template <> void FE<3,XYZ>::compute_constraints (DofConstraints &, DofMap &, const unsigned int, const Elem *) {}
944 
945 #endif // #ifdef LIBMESH_ENABLE_AMR
946 
947 // Full specialization of shapes_need_reinit() function for every dimension.
948 template <> bool FE<0,XYZ>::shapes_need_reinit() const { return true; }
949 template <> bool FE<1,XYZ>::shapes_need_reinit() const { return true; }
950 template <> bool FE<2,XYZ>::shapes_need_reinit() const { return true; }
951 template <> bool FE<3,XYZ>::shapes_need_reinit() const { return true; }
952 
953 
954 // Explicit instantiations for non-static FEXYZ member functions.
955 // These non-static member functions map more naturally to explicit
956 // instantiations than the functions above:
957 //
958 // 1.) Since they are member functions, they rely on
959 // private/protected member data, and therefore do not work well
960 // with the "anonymous function call" model we've used above for
961 // the specializations.
962 //
963 // 2.) There is (IMHO) less chance of the linker calling the
964 // wrong version of one of these member functions, since there is
965 // only one FEXYZ.
966 template void FEXYZ<0>::init_shape_functions(const std::vector<Point> &, const Elem *);
967 template void FEXYZ<1>::init_shape_functions(const std::vector<Point> &, const Elem *);
968 template void FEXYZ<2>::init_shape_functions(const std::vector<Point> &, const Elem *);
969 template void FEXYZ<3>::init_shape_functions(const std::vector<Point> &, const Elem *);
970 
971 template void FEXYZ<0>::compute_shape_functions(const Elem *,const std::vector<Point> &);
972 template void FEXYZ<1>::compute_shape_functions(const Elem *,const std::vector<Point> &);
973 template void FEXYZ<2>::compute_shape_functions(const Elem *,const std::vector<Point> &);
974 template void FEXYZ<3>::compute_shape_functions(const Elem *,const std::vector<Point> &);
975 
976 } // namespace libMesh
static unsigned int n_dofs(const ElemType t, const Order o)
virtual void init_shape_functions(const std::vector< Point > &qp, const Elem *e) libmesh_override
Update the various member data fields phi, dphidxi, dphideta, dphidzeta, etc.
Definition: fe_xyz.C:600
virtual bool shapes_need_reinit() const libmesh_override
unsigned int dim
ElemType
Defines an enum for geometric element types.
This is the base class from which all geometric element types are derived.
Definition: elem.h:89
static unsigned int n_dofs_at_node(const ElemType t, const Order o, const unsigned int n)
The libMesh namespace provides an interface to certain functionality in the library.
libmesh_assert(j)
This class handles the numbering of degrees of freedom on a mesh.
Definition: dof_map.h:167
A specific instantiation of the FEBase class.
Definition: fe.h:89
const dof_id_type n_nodes
Definition: tecplot_io.C:67
virtual FEContinuity get_continuity() const libmesh_override
static unsigned int n_shape_functions(const unsigned int dim, const FEType &fe_t, const ElemType t)
Definition: fe_interface.C:385
dof_id_type weight(const MeshBase &mesh, const processor_id_type pid)
Definition: mesh_tools.C:245
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:641
virtual bool is_hierarchic() const libmesh_override
PetscErrorCode Vec Mat libmesh_dbg_var(j)
static unsigned int n_dofs_per_elem(const ElemType t, const Order o)
FEContinuity
defines an enum for finite element types to libmesh_assert a certain level (or type? Hcurl?) of continuity.
std::string enum_to_string(const T e)
static void compute_constraints(DofConstraints &constraints, DofMap &dof_map, const unsigned int variable_number, const Elem *elem)
Computes the constraint matrix contributions (for non-conforming adapted meshes) corresponding to var...
static void nodal_soln(const Elem *elem, const Order o, const std::vector< Number > &elem_soln, std::vector< Number > &nodal_soln)
Build the nodal soln from the element soln.
Order
defines an enum for polynomial orders.
Definition: enum_order.h:32
virtual void compute_shape_functions(const Elem *elem, const std::vector< Point > &qp) libmesh_override
After having updated the jacobian and the transformation from local to global coordinates in FEAbstra...
Definition: fe_xyz.C:711
The constraint matrix storage format.
Definition: dof_map.h:96
static OutputShape shape_second_deriv(const ElemType t, const Order o, const unsigned int i, const unsigned int j, const Point &p)