fe_xyz.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2018 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 #include "libmesh/int_range.h"
27 
28 namespace libMesh
29 {
30 
31 // ------------------------------------------------------------
32 // XYZ-specific implementations
33 
34 // Anonymous namespace for local helper functions
35 namespace {
36 
37 void xyz_nodal_soln(const Elem * elem,
38  const Order order,
39  const std::vector<Number> & elem_soln,
40  std::vector<Number> & nodal_soln,
41  unsigned Dim)
42 {
43  const unsigned int n_nodes = elem->n_nodes();
44 
45  const ElemType elem_type = elem->type();
46 
47  nodal_soln.resize(n_nodes);
48 
49  const Order totalorder = static_cast<Order>(order + elem->p_level());
50 
51  switch (totalorder)
52  {
53  // Constant shape functions
54  case CONSTANT:
55  {
56  libmesh_assert_equal_to (elem_soln.size(), 1);
57 
58  const Number val = elem_soln[0];
59 
60  for (unsigned int n=0; n<n_nodes; n++)
61  nodal_soln[n] = val;
62 
63  return;
64  }
65 
66 
67  // For other orders do interpolation at the nodes
68  // explicitly.
69  default:
70  {
71  // FEType object to be passed to various FEInterface functions below.
72  FEType fe_type(totalorder, XYZ);
73 
74  const unsigned int n_sf =
75  // FE<Dim,T>::n_shape_functions(elem_type, totalorder);
76  FEInterface::n_shape_functions(Dim, fe_type, elem_type);
77 
78  for (unsigned int n=0; n<n_nodes; n++)
79  {
80  libmesh_assert_equal_to (elem_soln.size(), n_sf);
81 
82  // Zero before summation
83  nodal_soln[n] = 0;
84 
85  // u_i = Sum (alpha_i phi_i)
86  for (unsigned int i=0; i<n_sf; i++)
87  nodal_soln[n] += elem_soln[i] *
88  // FE<Dim,T>::shape(elem, order, i, elem->point(n));
89  FEInterface::shape(Dim, fe_type, elem, i, elem->point(n));
90  }
91 
92  return;
93  } // default
94  } // switch
95 } // xyz_nodal_soln()
96 
97 
98 
99 
100 
101 unsigned int xyz_n_dofs(const ElemType t, const Order o)
102 {
103  switch (o)
104  {
105 
106  // constant shape functions
107  // no matter what shape there is only one DOF.
108  case CONSTANT:
109  return (t != INVALID_ELEM) ? 1 : 0;
110 
111 
112  // Discontinuous linear shape functions
113  // expressed in the XYZ monomials.
114  case FIRST:
115  {
116  switch (t)
117  {
118  case NODEELEM:
119  return 1;
120 
121  case EDGE2:
122  case EDGE3:
123  case EDGE4:
124  return 2;
125 
126  case TRI3:
127  case TRISHELL3:
128  case TRI6:
129  case QUAD4:
130  case QUADSHELL4:
131  case QUAD8:
132  case QUADSHELL8:
133  case QUAD9:
134  return 3;
135 
136  case TET4:
137  case TET10:
138  case HEX8:
139  case HEX20:
140  case HEX27:
141  case PRISM6:
142  case PRISM15:
143  case PRISM18:
144  case PYRAMID5:
145  case PYRAMID13:
146  case PYRAMID14:
147  return 4;
148 
149  case INVALID_ELEM:
150  return 0;
151 
152  default:
153  libmesh_error_msg("ERROR: Bad ElemType = " << Utility::enum_to_string(t) << " for " << Utility::enum_to_string(o) << " order approximation!");
154  }
155  }
156 
157 
158  // Discontinuous quadratic shape functions
159  // expressed in the XYZ monomials.
160  case SECOND:
161  {
162  switch (t)
163  {
164  case NODEELEM:
165  return 1;
166 
167  case EDGE2:
168  case EDGE3:
169  case EDGE4:
170  return 3;
171 
172  case TRI3:
173  case TRISHELL3:
174  case TRI6:
175  case QUAD4:
176  case QUADSHELL4:
177  case QUAD8:
178  case QUADSHELL8:
179  case QUAD9:
180  return 6;
181 
182  case TET4:
183  case TET10:
184  case HEX8:
185  case HEX20:
186  case HEX27:
187  case PRISM6:
188  case PRISM15:
189  case PRISM18:
190  case PYRAMID5:
191  case PYRAMID13:
192  case PYRAMID14:
193  return 10;
194 
195  case INVALID_ELEM:
196  return 0;
197 
198  default:
199  libmesh_error_msg("ERROR: Bad ElemType = " << Utility::enum_to_string(t) << " for " << Utility::enum_to_string(o) << " order approximation!");
200  }
201  }
202 
203 
204  // Discontinuous cubic shape functions
205  // expressed in the XYZ monomials.
206  case THIRD:
207  {
208  switch (t)
209  {
210  case NODEELEM:
211  return 1;
212 
213  case EDGE2:
214  case EDGE3:
215  case EDGE4:
216  return 4;
217 
218  case TRI3:
219  case TRISHELL3:
220  case TRI6:
221  case QUAD4:
222  case QUADSHELL4:
223  case QUAD8:
224  case QUADSHELL8:
225  case QUAD9:
226  return 10;
227 
228  case TET4:
229  case TET10:
230  case HEX8:
231  case HEX20:
232  case HEX27:
233  case PRISM6:
234  case PRISM15:
235  case PRISM18:
236  case PYRAMID5:
237  case PYRAMID13:
238  case PYRAMID14:
239  return 20;
240 
241  case INVALID_ELEM:
242  return 0;
243 
244  default:
245  libmesh_error_msg("ERROR: Bad ElemType = " << Utility::enum_to_string(t) << " for " << Utility::enum_to_string(o) << " order approximation!");
246  }
247  }
248 
249 
250  // Discontinuous quartic shape functions
251  // expressed in the XYZ monomials.
252  case FOURTH:
253  {
254  switch (t)
255  {
256  case NODEELEM:
257  return 1;
258 
259  case EDGE2:
260  case EDGE3:
261  return 5;
262 
263  case TRI3:
264  case TRISHELL3:
265  case TRI6:
266  case QUAD4:
267  case QUADSHELL4:
268  case QUAD8:
269  case QUADSHELL8:
270  case QUAD9:
271  return 15;
272 
273  case TET4:
274  case TET10:
275  case HEX8:
276  case HEX20:
277  case HEX27:
278  case PRISM6:
279  case PRISM15:
280  case PRISM18:
281  case PYRAMID5:
282  case PYRAMID13:
283  case PYRAMID14:
284  return 35;
285 
286  case INVALID_ELEM:
287  return 0;
288 
289  default:
290  libmesh_error_msg("ERROR: Bad ElemType = " << Utility::enum_to_string(t) << " for " << Utility::enum_to_string(o) << " order approximation!");
291  }
292  }
293 
294 
295  default:
296  {
297  const unsigned int order = static_cast<unsigned int>(o);
298  switch (t)
299  {
300  case NODEELEM:
301  return 1;
302 
303  case EDGE2:
304  case EDGE3:
305  return (order+1);
306 
307  case TRI3:
308  case TRISHELL3:
309  case TRI6:
310  case QUAD4:
311  case QUADSHELL4:
312  case QUAD8:
313  case QUADSHELL8:
314  case QUAD9:
315  return (order+1)*(order+2)/2;
316 
317  case TET4:
318  case TET10:
319  case HEX8:
320  case HEX20:
321  case HEX27:
322  case PRISM6:
323  case PRISM15:
324  case PRISM18:
325  case PYRAMID5:
326  case PYRAMID13:
327  case PYRAMID14:
328  return (order+1)*(order+2)*(order+3)/6;
329 
330  case INVALID_ELEM:
331  return 0;
332 
333  default:
334  libmesh_error_msg("ERROR: Bad ElemType = " << Utility::enum_to_string(t) << " for " << Utility::enum_to_string(o) << " order approximation!");
335  }
336  }
337  }
338 }
339 
340 
341 
342 
343 unsigned int xyz_n_dofs_per_elem(const ElemType t,
344  const Order o)
345 {
346  switch (o)
347  {
348  // constant shape functions always have 1 DOF per element
349  case CONSTANT:
350  return (t != INVALID_ELEM) ? 1 : 0;
351 
352 
353  // Discontinuous linear shape functions
354  // expressed in the XYZ monomials.
355  case FIRST:
356  {
357  switch (t)
358  {
359  case NODEELEM:
360  return 1;
361 
362  // 1D linears have 2 DOFs per element
363  case EDGE2:
364  case EDGE3:
365  case EDGE4:
366  return 2;
367 
368  // 2D linears have 3 DOFs per element
369  case TRI3:
370  case TRISHELL3:
371  case TRI6:
372  case QUAD4:
373  case QUADSHELL4:
374  case QUAD8:
375  case QUADSHELL8:
376  case QUAD9:
377  return 3;
378 
379  // 3D linears have 4 DOFs per element
380  case TET4:
381  case TET10:
382  case HEX8:
383  case HEX20:
384  case HEX27:
385  case PRISM6:
386  case PRISM15:
387  case PRISM18:
388  case PYRAMID5:
389  case PYRAMID13:
390  case PYRAMID14:
391  return 4;
392 
393  case INVALID_ELEM:
394  return 0;
395 
396  default:
397  libmesh_error_msg("ERROR: Bad ElemType = " << Utility::enum_to_string(t) << " for " << Utility::enum_to_string(o) << " order approximation!");
398  }
399  }
400 
401 
402  // Discontinuous quadratic shape functions
403  // expressed in the XYZ monomials.
404  case SECOND:
405  {
406  switch (t)
407  {
408  case NODEELEM:
409  return 1;
410 
411  // 1D quadratics have 3 DOFs per element
412  case EDGE2:
413  case EDGE3:
414  case EDGE4:
415  return 3;
416 
417  // 2D quadratics have 6 DOFs per element
418  case TRI3:
419  case TRISHELL3:
420  case TRI6:
421  case QUAD4:
422  case QUADSHELL4:
423  case QUAD8:
424  case QUADSHELL8:
425  case QUAD9:
426  return 6;
427 
428  // 3D quadratics have 10 DOFs per element
429  case TET4:
430  case TET10:
431  case HEX8:
432  case HEX20:
433  case HEX27:
434  case PRISM6:
435  case PRISM15:
436  case PRISM18:
437  case PYRAMID5:
438  case PYRAMID13:
439  case PYRAMID14:
440  return 10;
441 
442  case INVALID_ELEM:
443  return 0;
444 
445  default:
446  libmesh_error_msg("ERROR: Bad ElemType = " << Utility::enum_to_string(t) << " for " << Utility::enum_to_string(o) << " order approximation!");
447  }
448  }
449 
450 
451  // Discontinuous cubic shape functions
452  // expressed in the XYZ monomials.
453  case THIRD:
454  {
455  switch (t)
456  {
457  case NODEELEM:
458  return 1;
459 
460  case EDGE2:
461  case EDGE3:
462  case EDGE4:
463  return 4;
464 
465  case TRI3:
466  case TRISHELL3:
467  case TRI6:
468  case QUAD4:
469  case QUADSHELL4:
470  case QUAD8:
471  case QUADSHELL8:
472  case QUAD9:
473  return 10;
474 
475  case TET4:
476  case TET10:
477  case HEX8:
478  case HEX20:
479  case HEX27:
480  case PRISM6:
481  case PRISM15:
482  case PRISM18:
483  case PYRAMID5:
484  case PYRAMID13:
485  case PYRAMID14:
486  return 20;
487 
488  case INVALID_ELEM:
489  return 0;
490 
491  default:
492  libmesh_error_msg("ERROR: Bad ElemType = " << Utility::enum_to_string(t) << " for " << Utility::enum_to_string(o) << " order approximation!");
493  }
494  }
495 
496 
497  // Discontinuous quartic shape functions
498  // expressed in the XYZ monomials.
499  case FOURTH:
500  {
501  switch (t)
502  {
503  case NODEELEM:
504  return 1;
505 
506  case EDGE2:
507  case EDGE3:
508  case EDGE4:
509  return 5;
510 
511  case TRI3:
512  case TRISHELL3:
513  case TRI6:
514  case QUAD4:
515  case QUADSHELL4:
516  case QUAD8:
517  case QUADSHELL8:
518  case QUAD9:
519  return 15;
520 
521  case TET4:
522  case TET10:
523  case HEX8:
524  case HEX20:
525  case HEX27:
526  case PRISM6:
527  case PRISM15:
528  case PRISM18:
529  case PYRAMID5:
530  case PYRAMID13:
531  case PYRAMID14:
532  return 35;
533 
534  case INVALID_ELEM:
535  return 0;
536 
537  default:
538  libmesh_error_msg("ERROR: Bad ElemType = " << Utility::enum_to_string(t) << " for " << Utility::enum_to_string(o) << " order approximation!");
539  }
540  }
541 
542  default:
543  {
544  const unsigned int order = static_cast<unsigned int>(o);
545  switch (t)
546  {
547  case NODEELEM:
548  return 1;
549 
550  case EDGE2:
551  case EDGE3:
552  return (order+1);
553 
554  case TRI3:
555  case TRISHELL3:
556  case TRI6:
557  case QUAD4:
558  case QUADSHELL4:
559  case QUAD8:
560  case QUADSHELL8:
561  case QUAD9:
562  return (order+1)*(order+2)/2;
563 
564  case TET4:
565  case TET10:
566  case HEX8:
567  case HEX20:
568  case HEX27:
569  case PRISM6:
570  case PRISM15:
571  case PRISM18:
572  case PYRAMID5:
573  case PYRAMID13:
574  case PYRAMID14:
575  return (order+1)*(order+2)*(order+3)/6;
576 
577  case INVALID_ELEM:
578  return 0;
579 
580  default:
581  libmesh_error_msg("ERROR: Bad ElemType = " << Utility::enum_to_string(t) << " for " << Utility::enum_to_string(o) << " order approximation!");
582  }
583  }
584  return 0;
585  }
586 }
587 
588 
589 } // anonymous namespace
590 
591 
592 
593 
594 
595 
596 
597 template <unsigned int Dim>
598 void FEXYZ<Dim>::init_shape_functions(const std::vector<Point> & qp,
599  const Elem * libmesh_dbg_var(elem))
600 {
601  libmesh_assert(elem);
602  this->calculations_started = true;
603 
604  // If the user forgot to request anything, we'll be safe and
605  // calculate everything:
606 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
607  if (!this->calculate_phi && !this->calculate_dphi && !this->calculate_d2phi)
608  this->calculate_phi = this->calculate_dphi = this->calculate_d2phi = true;
609 #else
610  if (!this->calculate_phi && !this->calculate_dphi)
611  this->calculate_phi = this->calculate_dphi = true;
612 #endif // LIBMESH_ENABLE_SECOND_DERIVATIVES
613 
614  // Start logging the shape function initialization
615  LOG_SCOPE("init_shape_functions()", "FE");
616 
617  // The number of quadrature points.
618  const std::size_t n_qp = qp.size();
619 
620  // Number of shape functions in the finite element approximation
621  // space.
622  const unsigned int n_approx_shape_functions =
623  this->n_shape_functions(this->get_type(),
624  this->get_order());
625 
626  // resize the vectors to hold current data
627  // Phi are the shape functions used for the FE approximation
628  {
629  // (note: GCC 3.4.0 requires the use of this-> here)
630  if (this->calculate_phi)
631  this->phi.resize (n_approx_shape_functions);
632  if (this->calculate_dphi)
633  {
634  this->dphi.resize (n_approx_shape_functions);
635  this->dphidx.resize (n_approx_shape_functions);
636  this->dphidy.resize (n_approx_shape_functions);
637  this->dphidz.resize (n_approx_shape_functions);
638  }
639 
640 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
641  if (this->calculate_d2phi)
642  {
643  this->d2phi.resize (n_approx_shape_functions);
644  this->d2phidx2.resize (n_approx_shape_functions);
645  this->d2phidxdy.resize (n_approx_shape_functions);
646  this->d2phidxdz.resize (n_approx_shape_functions);
647  this->d2phidy2.resize (n_approx_shape_functions);
648  this->d2phidydz.resize (n_approx_shape_functions);
649  this->d2phidz2.resize (n_approx_shape_functions);
650  this->d2phidxi2.resize (n_approx_shape_functions);
651  }
652 #endif // ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
653 
654  for (unsigned int i=0; i<n_approx_shape_functions; i++)
655  {
656  if (this->calculate_phi)
657  this->phi[i].resize (n_qp);
658  if (this->calculate_dphi)
659  {
660  this->dphi[i].resize (n_qp);
661  this->dphidx[i].resize (n_qp);
662  this->dphidy[i].resize (n_qp);
663  this->dphidz[i].resize (n_qp);
664  }
665 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
666  if (this->calculate_d2phi)
667  {
668  this->d2phi[i].resize (n_qp);
669  this->d2phidx2[i].resize (n_qp);
670  this->d2phidxdy[i].resize (n_qp);
671  this->d2phidxdz[i].resize (n_qp);
672  this->d2phidy2[i].resize (n_qp);
673  this->d2phidydz[i].resize (n_qp);
674  this->d2phidz2[i].resize (n_qp);
675  }
676 #endif // ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
677  }
678  }
679 
680 
681 
682 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
683  //------------------------------------------------------------
684  // Initialize the data fields, which should only be used for infinite
685  // elements, to some sensible values, so that using a FE with the
686  // variational formulation of an InfFE, correct element matrices are
687  // returned
688 
689  {
690  this->weight.resize (n_qp);
691  this->dweight.resize (n_qp);
692  this->dphase.resize (n_qp);
693 
694  for (unsigned int p=0; p<n_qp; p++)
695  {
696  this->weight[p] = 1.;
697  this->dweight[p].zero();
698  this->dphase[p].zero();
699  }
700 
701  }
702 #endif // ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
703 }
704 
705 
706 
707 
708 template <unsigned int Dim>
709 void FEXYZ<Dim>::compute_shape_functions (const Elem * elem,
710  const std::vector<Point> &)
711 {
712  libmesh_assert(elem);
713 
714  //-------------------------------------------------------------------------
715  // Compute the shape function values (and derivatives)
716  // at the Quadrature points. Note that the actual values
717  // have already been computed via init_shape_functions
718 
719  // Start logging the shape function computation
720  LOG_SCOPE("compute_shape_functions()", "FE");
721 
722  const std::vector<Point> & xyz_qp = this->get_xyz();
723 
724  // Compute the value of the derivative shape function i at quadrature point p
725  switch (this->dim)
726  {
727 
728  case 1:
729  {
730  if (this->calculate_phi)
731  for (auto i : index_range(this->phi))
732  for (auto p : index_range(this->phi[i]))
733  this->phi[i][p] = FE<Dim,XYZ>::shape (elem, this->fe_type.order, i, xyz_qp[p]);
734 
735  if (this->calculate_dphi)
736  for (auto i : index_range(this->dphi))
737  for (auto p : index_range(this->dphi[i]))
738  {
739  this->dphi[i][p](0) =
740  this->dphidx[i][p] = FE<Dim,XYZ>::shape_deriv (elem, this->fe_type.order, i, 0, xyz_qp[p]);
741 
742  this->dphi[i][p](1) = this->dphidy[i][p] = 0.;
743  this->dphi[i][p](2) = this->dphidz[i][p] = 0.;
744  }
745 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
746  if (this->calculate_d2phi)
747  for (auto i : index_range(this->d2phi))
748  for (auto p : index_range(this->d2phi[i]))
749  {
750  this->d2phi[i][p](0,0) =
751  this->d2phidx2[i][p] = FE<Dim,XYZ>::shape_second_deriv (elem, this->fe_type.order, i, 0, xyz_qp[p]);
752 
753 #if LIBMESH_DIM>1
754  this->d2phi[i][p](0,1) = this->d2phidxdy[i][p] =
755  this->d2phi[i][p](1,0) = 0.;
756  this->d2phi[i][p](1,1) = this->d2phidy2[i][p] = 0.;
757 #if LIBMESH_DIM>2
758  this->d2phi[i][p](0,2) = this->d2phidxdz[i][p] =
759  this->d2phi[i][p](2,0) = 0.;
760  this->d2phi[i][p](1,2) = this->d2phidydz[i][p] =
761  this->d2phi[i][p](2,1) = 0.;
762  this->d2phi[i][p](2,2) = this->d2phidz2[i][p] = 0.;
763 #endif
764 #endif
765  }
766 #endif // ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
767 
768  // All done
769  break;
770  }
771 
772  case 2:
773  {
774  if (this->calculate_phi)
775  for (auto i : index_range(this->phi))
776  for (auto p : index_range(this->phi[i]))
777  this->phi[i][p] = FE<Dim,XYZ>::shape (elem, this->fe_type.order, i, xyz_qp[p]);
778 
779  if (this->calculate_dphi)
780  for (auto i : index_range(this->dphi))
781  for (auto p : index_range(this->dphi[i]))
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 (auto i : index_range(this->d2phi))
797  for (auto p : index_range(this->d2phi[i]))
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 (auto i : index_range(this->phi))
824  for (auto p : index_range(this->phi[i]))
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 (auto i : index_range(this->dphi))
829  for (auto p : index_range(this->dphi[i]))
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 (auto i : index_range(this->d2phi))
843  for (auto p : index_range(this->d2phi[i]))
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
libMesh::FE::n_dofs
static unsigned int n_dofs(const ElemType t, const Order o)
libMesh::FE::shape
static OutputShape shape(const ElemType t, const Order o, const unsigned int i, const Point &p)
libMesh::Elem
The base class for all geometric element types.
Definition: elem.h:100
libMesh::index_range
IntRange< std::size_t > index_range(const std::vector< T > &vec)
Definition: int_range.h:104
libMesh::FE::n_dofs_at_node
static unsigned int n_dofs_at_node(const ElemType t, const Order o, const unsigned int n)
fe_interface.h
libMesh::FE::shapes_need_reinit
virtual bool shapes_need_reinit() const override
libMesh::FEXYZ::init_shape_functions
virtual void init_shape_functions(const std::vector< Point > &qp, const Elem *e) override
Definition: fe_xyz.C:598
libMesh::FE::is_hierarchic
virtual bool is_hierarchic() const override
libMesh::DofMap
Manages the degrees of freedom (DOFs) in a simulation.
Definition: dof_map.h:176
libMesh::FE
Template class which generates the different FE families and orders.
Definition: fe.h:89
n_nodes
const dof_id_type n_nodes
Definition: tecplot_io.C:68
libMesh::FEInterface::n_shape_functions
static unsigned int n_shape_functions(const unsigned int dim, const FEType &fe_t, const ElemType t)
Definition: fe_interface.C:428
libMesh::MeshTools::weight
dof_id_type weight(const MeshBase &mesh, const processor_id_type pid)
Definition: mesh_tools.C:233
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:657
libmesh_logging.h
int_range.h
libMesh::FE::n_dofs_per_elem
static unsigned int n_dofs_per_elem(const ElemType t, const Order o)
libMesh::FE::get_continuity
virtual FEContinuity get_continuity() const override
libMesh::Utility::enum_to_string
std::string enum_to_string(const T e)
libMesh::FE::compute_constraints
static void compute_constraints(DofConstraints &constraints, DofMap &dof_map, const unsigned int variable_number, const Elem *elem)
string_to_enum.h
libMesh::FEXYZ::compute_shape_functions
virtual void compute_shape_functions(const Elem *elem, const std::vector< Point > &qp) override
Definition: fe_xyz.C:709
libMesh::FE::nodal_soln
static void nodal_soln(const Elem *elem, const Order o, const std::vector< Number > &elem_soln, std::vector< Number > &nodal_soln)
libMesh::FE::shape_second_deriv
static OutputShape shape_second_deriv(const ElemType t, const Order o, const unsigned int i, const unsigned int j, const Point &p)
generated by doxygen