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

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

Hosted By:
SourceForge.net Logo