fe_abstract.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/fe.h"
22 #include "libmesh/libmesh_logging.h"
23 #include "libmesh/enum_elem_type.h"
24 
25 // For projection code:
26 #include "libmesh/boundary_info.h"
27 #include "libmesh/mesh_base.h"
28 #include "libmesh/dense_matrix.h"
29 #include "libmesh/dense_vector.h"
30 #include "libmesh/dof_map.h"
31 #include "libmesh/elem.h"
32 #include "libmesh/fe_interface.h"
33 #include "libmesh/numeric_vector.h"
34 #include "libmesh/periodic_boundaries.h"
35 #include "libmesh/periodic_boundary.h"
36 #include "libmesh/quadrature.h"
37 #include "libmesh/quadrature_gauss.h"
38 #include "libmesh/remote_elem.h"
39 #include "libmesh/tensor_value.h"
40 #include "libmesh/threads.h"
41 #include "libmesh/enum_elem_type.h"
42 
43 namespace libMesh
44 {
45 
46 FEAbstract::FEAbstract(const unsigned int d,
47  const FEType & fet) :
48  _fe_map( FEMap::build(fet) ),
49  dim(d),
50  calculations_started(false),
51  calculate_phi(false),
52  calculate_dphi(false),
53  calculate_d2phi(false),
54  calculate_curl_phi(false),
55  calculate_div_phi(false),
56  calculate_dphiref(false),
57  fe_type(fet),
58  elem_type(INVALID_ELEM),
59  _p_level(0),
60  qrule(nullptr),
61  shapes_on_quadrature(false)
62 {
63 }
64 
65 
66 FEAbstract::~FEAbstract()
67 {
68 }
69 
70 
71 std::unique_ptr<FEAbstract> FEAbstract::build(const unsigned int dim,
72  const FEType & fet)
73 {
74  switch (dim)
75  {
76  // 0D
77  case 0:
78  {
79  switch (fet.family)
80  {
81  case CLOUGH:
82  return libmesh_make_unique<FE<0,CLOUGH>>(fet);
83 
84  case HERMITE:
85  return libmesh_make_unique<FE<0,HERMITE>>(fet);
86 
87  case LAGRANGE:
88  return libmesh_make_unique<FE<0,LAGRANGE>>(fet);
89 
90  case LAGRANGE_VEC:
91  return libmesh_make_unique<FE<0,LAGRANGE_VEC>>(fet);
92 
93  case L2_LAGRANGE:
94  return libmesh_make_unique<FE<0,L2_LAGRANGE>>(fet);
95 
96  case HIERARCHIC:
97  return libmesh_make_unique<FE<0,HIERARCHIC>>(fet);
98 
99  case L2_HIERARCHIC:
100  return libmesh_make_unique<FE<0,L2_HIERARCHIC>>(fet);
101 
102  case MONOMIAL:
103  return libmesh_make_unique<FE<0,MONOMIAL>>(fet);
104 
105 #ifdef LIBMESH_ENABLE_HIGHER_ORDER_SHAPES
106  case SZABAB:
107  return libmesh_make_unique<FE<0,SZABAB>>(fet);
108 
109  case BERNSTEIN:
110  return libmesh_make_unique<FE<0,BERNSTEIN>>(fet);
111 #endif
112 
113  case XYZ:
114  return libmesh_make_unique<FEXYZ<0>>(fet);
115 
116  case SCALAR:
117  return libmesh_make_unique<FEScalar<0>>(fet);
118 
119  default:
120  libmesh_error_msg("ERROR: Bad FEType.family= " << fet.family);
121  }
122  }
123  // 1D
124  case 1:
125  {
126  switch (fet.family)
127  {
128  case CLOUGH:
129  return libmesh_make_unique<FE<1,CLOUGH>>(fet);
130 
131  case HERMITE:
132  return libmesh_make_unique<FE<1,HERMITE>>(fet);
133 
134  case LAGRANGE:
135  return libmesh_make_unique<FE<1,LAGRANGE>>(fet);
136 
137  case LAGRANGE_VEC:
138  return libmesh_make_unique<FE<1,LAGRANGE_VEC>>(fet);
139 
140  case L2_LAGRANGE:
141  return libmesh_make_unique<FE<1,L2_LAGRANGE>>(fet);
142 
143  case HIERARCHIC:
144  return libmesh_make_unique<FE<1,HIERARCHIC>>(fet);
145 
146  case L2_HIERARCHIC:
147  return libmesh_make_unique<FE<1,L2_HIERARCHIC>>(fet);
148 
149  case MONOMIAL:
150  return libmesh_make_unique<FE<1,MONOMIAL>>(fet);
151 
152 #ifdef LIBMESH_ENABLE_HIGHER_ORDER_SHAPES
153  case SZABAB:
154  return libmesh_make_unique<FE<1,SZABAB>>(fet);
155 
156  case BERNSTEIN:
157  return libmesh_make_unique<FE<1,BERNSTEIN>>(fet);
158 #endif
159 
160  case XYZ:
161  return libmesh_make_unique<FEXYZ<1>>(fet);
162 
163  case SCALAR:
164  return libmesh_make_unique<FEScalar<1>>(fet);
165 
166  default:
167  libmesh_error_msg("ERROR: Bad FEType.family= " << fet.family);
168  }
169  }
170 
171 
172  // 2D
173  case 2:
174  {
175  switch (fet.family)
176  {
177  case CLOUGH:
178  return libmesh_make_unique<FE<2,CLOUGH>>(fet);
179 
180  case HERMITE:
181  return libmesh_make_unique<FE<2,HERMITE>>(fet);
182 
183  case LAGRANGE:
184  return libmesh_make_unique<FE<2,LAGRANGE>>(fet);
185 
186  case LAGRANGE_VEC:
187  return libmesh_make_unique<FE<2,LAGRANGE_VEC>>(fet);
188 
189  case L2_LAGRANGE:
190  return libmesh_make_unique<FE<2,L2_LAGRANGE>>(fet);
191 
192  case HIERARCHIC:
193  return libmesh_make_unique<FE<2,HIERARCHIC>>(fet);
194 
195  case L2_HIERARCHIC:
196  return libmesh_make_unique<FE<2,L2_HIERARCHIC>>(fet);
197 
198  case MONOMIAL:
199  return libmesh_make_unique<FE<2,MONOMIAL>>(fet);
200 
201 #ifdef LIBMESH_ENABLE_HIGHER_ORDER_SHAPES
202  case SZABAB:
203  return libmesh_make_unique<FE<2,SZABAB>>(fet);
204 
205  case BERNSTEIN:
206  return libmesh_make_unique<FE<2,BERNSTEIN>>(fet);
207 #endif
208 
209  case XYZ:
210  return libmesh_make_unique<FEXYZ<2>>(fet);
211 
212  case SCALAR:
213  return libmesh_make_unique<FEScalar<2>>(fet);
214 
215  case NEDELEC_ONE:
216  return libmesh_make_unique<FENedelecOne<2>>(fet);
217 
218  case SUBDIVISION:
219  return libmesh_make_unique<FESubdivision>(fet);
220 
221  default:
222  libmesh_error_msg("ERROR: Bad FEType.family= " << fet.family);
223  }
224  }
225 
226 
227  // 3D
228  case 3:
229  {
230  switch (fet.family)
231  {
232  case CLOUGH:
233  libmesh_error_msg("ERROR: Clough-Tocher elements currently only support 1D and 2D");
234 
235  case HERMITE:
236  return libmesh_make_unique<FE<3,HERMITE>>(fet);
237 
238  case LAGRANGE:
239  return libmesh_make_unique<FE<3,LAGRANGE>>(fet);
240 
241  case LAGRANGE_VEC:
242  return libmesh_make_unique<FE<3,LAGRANGE_VEC>>(fet);
243 
244  case L2_LAGRANGE:
245  return libmesh_make_unique<FE<3,L2_LAGRANGE>>(fet);
246 
247  case HIERARCHIC:
248  return libmesh_make_unique<FE<3,HIERARCHIC>>(fet);
249 
250  case L2_HIERARCHIC:
251  return libmesh_make_unique<FE<3,L2_HIERARCHIC>>(fet);
252 
253  case MONOMIAL:
254  return libmesh_make_unique<FE<3,MONOMIAL>>(fet);
255 
256 #ifdef LIBMESH_ENABLE_HIGHER_ORDER_SHAPES
257  case SZABAB:
258  return libmesh_make_unique<FE<3,SZABAB>>(fet);
259 
260  case BERNSTEIN:
261  return libmesh_make_unique<FE<3,BERNSTEIN>>(fet);
262 #endif
263 
264  case XYZ:
265  return libmesh_make_unique<FEXYZ<3>>(fet);
266 
267  case SCALAR:
268  return libmesh_make_unique<FEScalar<3>>(fet);
269 
270  case NEDELEC_ONE:
271  return libmesh_make_unique<FENedelecOne<3>>(fet);
272 
273  default:
274  libmesh_error_msg("ERROR: Bad FEType.family= " << fet.family);
275  }
276  }
277 
278  default:
279  libmesh_error_msg("Invalid dimension dim = " << dim);
280  }
281 }
282 
283 void FEAbstract::get_refspace_nodes(const ElemType itemType, std::vector<Point> & nodes)
284 {
285  switch(itemType)
286  {
287  case EDGE2:
288  {
289  nodes.resize(2);
290  nodes[0] = Point (-1.,0.,0.);
291  nodes[1] = Point (1.,0.,0.);
292  return;
293  }
294  case EDGE3:
295  {
296  nodes.resize(3);
297  nodes[0] = Point (-1.,0.,0.);
298  nodes[1] = Point (1.,0.,0.);
299  nodes[2] = Point (0.,0.,0.);
300  return;
301  }
302  case TRI3:
303  case TRISHELL3:
304  {
305  nodes.resize(3);
306  nodes[0] = Point (0.,0.,0.);
307  nodes[1] = Point (1.,0.,0.);
308  nodes[2] = Point (0.,1.,0.);
309  return;
310  }
311  case TRI6:
312  {
313  nodes.resize(6);
314  nodes[0] = Point (0.,0.,0.);
315  nodes[1] = Point (1.,0.,0.);
316  nodes[2] = Point (0.,1.,0.);
317  nodes[3] = Point (.5,0.,0.);
318  nodes[4] = Point (.5,.5,0.);
319  nodes[5] = Point (0.,.5,0.);
320  return;
321  }
322  case QUAD4:
323  case QUADSHELL4:
324  {
325  nodes.resize(4);
326  nodes[0] = Point (-1.,-1.,0.);
327  nodes[1] = Point (1.,-1.,0.);
328  nodes[2] = Point (1.,1.,0.);
329  nodes[3] = Point (-1.,1.,0.);
330  return;
331  }
332  case QUAD8:
333  case QUADSHELL8:
334  {
335  nodes.resize(8);
336  nodes[0] = Point (-1.,-1.,0.);
337  nodes[1] = Point (1.,-1.,0.);
338  nodes[2] = Point (1.,1.,0.);
339  nodes[3] = Point (-1.,1.,0.);
340  nodes[4] = Point (0.,-1.,0.);
341  nodes[5] = Point (1.,0.,0.);
342  nodes[6] = Point (0.,1.,0.);
343  nodes[7] = Point (-1.,0.,0.);
344  return;
345  }
346  case QUAD9:
347  {
348  nodes.resize(9);
349  nodes[0] = Point (-1.,-1.,0.);
350  nodes[1] = Point (1.,-1.,0.);
351  nodes[2] = Point (1.,1.,0.);
352  nodes[3] = Point (-1.,1.,0.);
353  nodes[4] = Point (0.,-1.,0.);
354  nodes[5] = Point (1.,0.,0.);
355  nodes[6] = Point (0.,1.,0.);
356  nodes[7] = Point (-1.,0.,0.);
357  nodes[8] = Point (0.,0.,0.);
358  return;
359  }
360  case TET4:
361  {
362  nodes.resize(4);
363  nodes[0] = Point (0.,0.,0.);
364  nodes[1] = Point (1.,0.,0.);
365  nodes[2] = Point (0.,1.,0.);
366  nodes[3] = Point (0.,0.,1.);
367  return;
368  }
369  case TET10:
370  {
371  nodes.resize(10);
372  nodes[0] = Point (0.,0.,0.);
373  nodes[1] = Point (1.,0.,0.);
374  nodes[2] = Point (0.,1.,0.);
375  nodes[3] = Point (0.,0.,1.);
376  nodes[4] = Point (.5,0.,0.);
377  nodes[5] = Point (.5,.5,0.);
378  nodes[6] = Point (0.,.5,0.);
379  nodes[7] = Point (0.,0.,.5);
380  nodes[8] = Point (.5,0.,.5);
381  nodes[9] = Point (0.,.5,.5);
382  return;
383  }
384  case HEX8:
385  {
386  nodes.resize(8);
387  nodes[0] = Point (-1.,-1.,-1.);
388  nodes[1] = Point (1.,-1.,-1.);
389  nodes[2] = Point (1.,1.,-1.);
390  nodes[3] = Point (-1.,1.,-1.);
391  nodes[4] = Point (-1.,-1.,1.);
392  nodes[5] = Point (1.,-1.,1.);
393  nodes[6] = Point (1.,1.,1.);
394  nodes[7] = Point (-1.,1.,1.);
395  return;
396  }
397  case HEX20:
398  {
399  nodes.resize(20);
400  nodes[0] = Point (-1.,-1.,-1.);
401  nodes[1] = Point (1.,-1.,-1.);
402  nodes[2] = Point (1.,1.,-1.);
403  nodes[3] = Point (-1.,1.,-1.);
404  nodes[4] = Point (-1.,-1.,1.);
405  nodes[5] = Point (1.,-1.,1.);
406  nodes[6] = Point (1.,1.,1.);
407  nodes[7] = Point (-1.,1.,1.);
408  nodes[8] = Point (0.,-1.,-1.);
409  nodes[9] = Point (1.,0.,-1.);
410  nodes[10] = Point (0.,1.,-1.);
411  nodes[11] = Point (-1.,0.,-1.);
412  nodes[12] = Point (-1.,-1.,0.);
413  nodes[13] = Point (1.,-1.,0.);
414  nodes[14] = Point (1.,1.,0.);
415  nodes[15] = Point (-1.,1.,0.);
416  nodes[16] = Point (0.,-1.,1.);
417  nodes[17] = Point (1.,0.,1.);
418  nodes[18] = Point (0.,1.,1.);
419  nodes[19] = Point (-1.,0.,1.);
420  return;
421  }
422  case HEX27:
423  {
424  nodes.resize(27);
425  nodes[0] = Point (-1.,-1.,-1.);
426  nodes[1] = Point (1.,-1.,-1.);
427  nodes[2] = Point (1.,1.,-1.);
428  nodes[3] = Point (-1.,1.,-1.);
429  nodes[4] = Point (-1.,-1.,1.);
430  nodes[5] = Point (1.,-1.,1.);
431  nodes[6] = Point (1.,1.,1.);
432  nodes[7] = Point (-1.,1.,1.);
433  nodes[8] = Point (0.,-1.,-1.);
434  nodes[9] = Point (1.,0.,-1.);
435  nodes[10] = Point (0.,1.,-1.);
436  nodes[11] = Point (-1.,0.,-1.);
437  nodes[12] = Point (-1.,-1.,0.);
438  nodes[13] = Point (1.,-1.,0.);
439  nodes[14] = Point (1.,1.,0.);
440  nodes[15] = Point (-1.,1.,0.);
441  nodes[16] = Point (0.,-1.,1.);
442  nodes[17] = Point (1.,0.,1.);
443  nodes[18] = Point (0.,1.,1.);
444  nodes[19] = Point (-1.,0.,1.);
445  nodes[20] = Point (0.,0.,-1.);
446  nodes[21] = Point (0.,-1.,0.);
447  nodes[22] = Point (1.,0.,0.);
448  nodes[23] = Point (0.,1.,0.);
449  nodes[24] = Point (-1.,0.,0.);
450  nodes[25] = Point (0.,0.,1.);
451  nodes[26] = Point (0.,0.,0.);
452  return;
453  }
454  case PRISM6:
455  {
456  nodes.resize(6);
457  nodes[0] = Point (0.,0.,-1.);
458  nodes[1] = Point (1.,0.,-1.);
459  nodes[2] = Point (0.,1.,-1.);
460  nodes[3] = Point (0.,0.,1.);
461  nodes[4] = Point (1.,0.,1.);
462  nodes[5] = Point (0.,1.,1.);
463  return;
464  }
465  case PRISM15:
466  {
467  nodes.resize(15);
468  nodes[0] = Point (0.,0.,-1.);
469  nodes[1] = Point (1.,0.,-1.);
470  nodes[2] = Point (0.,1.,-1.);
471  nodes[3] = Point (0.,0.,1.);
472  nodes[4] = Point (1.,0.,1.);
473  nodes[5] = Point (0.,1.,1.);
474  nodes[6] = Point (.5,0.,-1.);
475  nodes[7] = Point (.5,.5,-1.);
476  nodes[8] = Point (0.,.5,-1.);
477  nodes[9] = Point (0.,0.,0.);
478  nodes[10] = Point (1.,0.,0.);
479  nodes[11] = Point (0.,1.,0.);
480  nodes[12] = Point (.5,0.,1.);
481  nodes[13] = Point (.5,.5,1.);
482  nodes[14] = Point (0.,.5,1.);
483  return;
484  }
485  case PRISM18:
486  {
487  nodes.resize(18);
488  nodes[0] = Point (0.,0.,-1.);
489  nodes[1] = Point (1.,0.,-1.);
490  nodes[2] = Point (0.,1.,-1.);
491  nodes[3] = Point (0.,0.,1.);
492  nodes[4] = Point (1.,0.,1.);
493  nodes[5] = Point (0.,1.,1.);
494  nodes[6] = Point (.5,0.,-1.);
495  nodes[7] = Point (.5,.5,-1.);
496  nodes[8] = Point (0.,.5,-1.);
497  nodes[9] = Point (0.,0.,0.);
498  nodes[10] = Point (1.,0.,0.);
499  nodes[11] = Point (0.,1.,0.);
500  nodes[12] = Point (.5,0.,1.);
501  nodes[13] = Point (.5,.5,1.);
502  nodes[14] = Point (0.,.5,1.);
503  nodes[15] = Point (.5,0.,0.);
504  nodes[16] = Point (.5,.5,0.);
505  nodes[17] = Point (0.,.5,0.);
506  return;
507  }
508  case PYRAMID5:
509  {
510  nodes.resize(5);
511  nodes[0] = Point (-1.,-1.,0.);
512  nodes[1] = Point (1.,-1.,0.);
513  nodes[2] = Point (1.,1.,0.);
514  nodes[3] = Point (-1.,1.,0.);
515  nodes[4] = Point (0.,0.,1.);
516  return;
517  }
518  case PYRAMID13:
519  {
520  nodes.resize(13);
521 
522  // base corners
523  nodes[0] = Point (-1.,-1.,0.);
524  nodes[1] = Point (1.,-1.,0.);
525  nodes[2] = Point (1.,1.,0.);
526  nodes[3] = Point (-1.,1.,0.);
527 
528  // apex
529  nodes[4] = Point (0.,0.,1.);
530 
531  // base midedge
532  nodes[5] = Point (0.,-1.,0.);
533  nodes[6] = Point (1.,0.,0.);
534  nodes[7] = Point (0.,1.,0.);
535  nodes[8] = Point (-1,0.,0.);
536 
537  // lateral midedge
538  nodes[9] = Point (-.5,-.5,.5);
539  nodes[10] = Point (.5,-.5,.5);
540  nodes[11] = Point (.5,.5,.5);
541  nodes[12] = Point (-.5,.5,.5);
542 
543  return;
544  }
545  case PYRAMID14:
546  {
547  nodes.resize(14);
548 
549  // base corners
550  nodes[0] = Point (-1.,-1.,0.);
551  nodes[1] = Point (1.,-1.,0.);
552  nodes[2] = Point (1.,1.,0.);
553  nodes[3] = Point (-1.,1.,0.);
554 
555  // apex
556  nodes[4] = Point (0.,0.,1.);
557 
558  // base midedge
559  nodes[5] = Point (0.,-1.,0.);
560  nodes[6] = Point (1.,0.,0.);
561  nodes[7] = Point (0.,1.,0.);
562  nodes[8] = Point (-1,0.,0.);
563 
564  // lateral midedge
565  nodes[9] = Point (-.5,-.5,.5);
566  nodes[10] = Point (.5,-.5,.5);
567  nodes[11] = Point (.5,.5,.5);
568  nodes[12] = Point (-.5,.5,.5);
569 
570  // base center
571  nodes[13] = Point (0.,0.,0.);
572 
573  return;
574  }
575 
576  default:
577  libmesh_error_msg("ERROR: Unknown element type " << itemType);
578  }
579 }
580 
581 bool FEAbstract::on_reference_element(const Point & p, const ElemType t, const Real eps)
582 {
583  libmesh_assert_greater_equal (eps, 0.);
584 
585  const Real xi = p(0);
586 #if LIBMESH_DIM > 1
587  const Real eta = p(1);
588 #else
589  const Real eta = 0.;
590 #endif
591 #if LIBMESH_DIM > 2
592  const Real zeta = p(2);
593 #else
594  const Real zeta = 0.;
595 #endif
596 
597  switch (t)
598  {
599  case NODEELEM:
600  {
601  return (!xi && !eta && !zeta);
602  }
603  case EDGE2:
604  case EDGE3:
605  case EDGE4:
606  {
607  // The reference 1D element is [-1,1].
608  if ((xi >= -1.-eps) &&
609  (xi <= 1.+eps))
610  return true;
611 
612  return false;
613  }
614 
615 
616  case TRI3:
617  case TRISHELL3:
618  case TRI6:
619  {
620  // The reference triangle is isosceles
621  // and is bound by xi=0, eta=0, and xi+eta=1.
622  if ((xi >= 0.-eps) &&
623  (eta >= 0.-eps) &&
624  ((xi + eta) <= 1.+eps))
625  return true;
626 
627  return false;
628  }
629 
630 
631  case QUAD4:
632  case QUADSHELL4:
633  case QUAD8:
634  case QUADSHELL8:
635  case QUAD9:
636  {
637  // The reference quadrilateral element is [-1,1]^2.
638  if ((xi >= -1.-eps) &&
639  (xi <= 1.+eps) &&
640  (eta >= -1.-eps) &&
641  (eta <= 1.+eps))
642  return true;
643 
644  return false;
645  }
646 
647 
648  case TET4:
649  case TET10:
650  {
651  // The reference tetrahedral is isosceles
652  // and is bound by xi=0, eta=0, zeta=0,
653  // and xi+eta+zeta=1.
654  if ((xi >= 0.-eps) &&
655  (eta >= 0.-eps) &&
656  (zeta >= 0.-eps) &&
657  ((xi + eta + zeta) <= 1.+eps))
658  return true;
659 
660  return false;
661  }
662 
663 
664  case HEX8:
665  case HEX20:
666  case HEX27:
667  {
668  /*
669  if ((xi >= -1.) &&
670  (xi <= 1.) &&
671  (eta >= -1.) &&
672  (eta <= 1.) &&
673  (zeta >= -1.) &&
674  (zeta <= 1.))
675  return true;
676  */
677 
678  // The reference hexahedral element is [-1,1]^3.
679  if ((xi >= -1.-eps) &&
680  (xi <= 1.+eps) &&
681  (eta >= -1.-eps) &&
682  (eta <= 1.+eps) &&
683  (zeta >= -1.-eps) &&
684  (zeta <= 1.+eps))
685  {
686  // libMesh::out << "Strange Point:\n";
687  // p.print();
688  return true;
689  }
690 
691  return false;
692  }
693 
694  case PRISM6:
695  case PRISM15:
696  case PRISM18:
697  {
698  // Figure this one out...
699  // inside the reference triangle with zeta in [-1,1]
700  if ((xi >= 0.-eps) &&
701  (eta >= 0.-eps) &&
702  (zeta >= -1.-eps) &&
703  (zeta <= 1.+eps) &&
704  ((xi + eta) <= 1.+eps))
705  return true;
706 
707  return false;
708  }
709 
710 
711  case PYRAMID5:
712  case PYRAMID13:
713  case PYRAMID14:
714  {
715  // Check that the point is on the same side of all the faces
716  // by testing whether:
717  //
718  // n_i.(x - x_i) <= 0
719  //
720  // for each i, where:
721  // n_i is the outward normal of face i,
722  // x_i is a point on face i.
723  if ((-eta - 1. + zeta <= 0.+eps) &&
724  ( xi - 1. + zeta <= 0.+eps) &&
725  ( eta - 1. + zeta <= 0.+eps) &&
726  ( -xi - 1. + zeta <= 0.+eps) &&
727  ( zeta >= 0.-eps))
728  return true;
729 
730  return false;
731  }
732 
733 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
734  case INFHEX8:
735  case INFHEX16:
736  case INFHEX18:
737  {
738  // The reference infhex8 is a [-1,1]^3.
739  if ((xi >= -1.-eps) &&
740  (xi <= 1.+eps) &&
741  (eta >= -1.-eps) &&
742  (eta <= 1.+eps) &&
743  (zeta >= -1.-eps) &&
744  (zeta <= 1.+eps))
745  {
746  return true;
747  }
748  return false;
749  }
750 
751  case INFPRISM6:
752  case INFPRISM12:
753  {
754  // inside the reference triangle with zeta in [-1,1]
755  if ((xi >= 0.-eps) &&
756  (eta >= 0.-eps) &&
757  (zeta >= -1.-eps) &&
758  (zeta <= 1.+eps) &&
759  ((xi + eta) <= 1.+eps))
760  {
761  return true;
762  }
763 
764  return false;
765  }
766 #endif
767 
768  default:
769  libmesh_error_msg("ERROR: Unknown element type " << t);
770  }
771 
772  // If we get here then the point is _not_ in the
773  // reference element. Better return false.
774 
775  return false;
776 }
777 
778 
779 
780 void FEAbstract::print_JxW(std::ostream & os) const
781 {
782  this->_fe_map->print_JxW(os);
783 }
784 
785 
786 
787 void FEAbstract::print_xyz(std::ostream & os) const
788 {
789  this->_fe_map->print_xyz(os);
790 }
791 
792 
793 void FEAbstract::print_info(std::ostream & os) const
794 {
795  os << "phi[i][j]: Shape function i at quadrature pt. j" << std::endl;
796  this->print_phi(os);
797 
798  os << "dphi[i][j]: Shape function i's gradient at quadrature pt. j" << std::endl;
799  this->print_dphi(os);
800 
801  os << "XYZ locations of the quadrature pts." << std::endl;
802  this->print_xyz(os);
803 
804  os << "Values of JxW at the quadrature pts." << std::endl;
805  this->print_JxW(os);
806 }
807 
808 
809 std::ostream & operator << (std::ostream & os, const FEAbstract & fe)
810 {
811  fe.print_info(os);
812  return os;
813 }
814 
815 
816 
817 #ifdef LIBMESH_ENABLE_AMR
818 
819 #ifdef LIBMESH_ENABLE_NODE_CONSTRAINTS
820 void FEAbstract::compute_node_constraints (NodeConstraints & constraints,
821  const Elem * elem)
822 {
823  libmesh_assert(elem);
824 
825  const unsigned int Dim = elem->dim();
826 
827  // Only constrain elements in 2,3D.
828  if (Dim == 1)
829  return;
830 
831  // Only constrain active and ancestor elements
832  if (elem->subactive())
833  return;
834 
835  // We currently always use LAGRANGE mappings for geometry
836  const FEType fe_type(elem->default_order(), LAGRANGE);
837 
838  // Pull objects out of the loop to reduce heap operations
839  std::vector<const Node *> my_nodes, parent_nodes;
840  std::unique_ptr<const Elem> my_side, parent_side;
841 
842  // Look at the element faces. Check to see if we need to
843  // build constraints.
844  for (auto s : elem->side_index_range())
845  if (elem->neighbor_ptr(s) != nullptr &&
846  elem->neighbor_ptr(s) != remote_elem)
847  if (elem->neighbor_ptr(s)->level() < elem->level()) // constrain dofs shared between
848  { // this element and ones coarser
849  // than this element.
850  // Get pointers to the elements of interest and its parent.
851  const Elem * parent = elem->parent();
852 
853  // This can't happen... Only level-0 elements have nullptr
854  // parents, and no level-0 elements can be at a higher
855  // level than their neighbors!
856  libmesh_assert(parent);
857 
858  elem->build_side_ptr(my_side, s);
859  parent->build_side_ptr(parent_side, s);
860 
861  const unsigned int n_side_nodes = my_side->n_nodes();
862 
863  my_nodes.clear();
864  my_nodes.reserve (n_side_nodes);
865  parent_nodes.clear();
866  parent_nodes.reserve (n_side_nodes);
867 
868  for (unsigned int n=0; n != n_side_nodes; ++n)
869  my_nodes.push_back(my_side->node_ptr(n));
870 
871  for (unsigned int n=0; n != n_side_nodes; ++n)
872  parent_nodes.push_back(parent_side->node_ptr(n));
873 
874  for (unsigned int my_side_n=0;
875  my_side_n < n_side_nodes;
876  my_side_n++)
877  {
878  libmesh_assert_less (my_side_n, FEInterface::n_dofs(Dim-1, fe_type, my_side->type()));
879 
880  const Node * my_node = my_nodes[my_side_n];
881 
882  // The support point of the DOF
883  const Point & support_point = *my_node;
884 
885  // Figure out where my node lies on their reference element.
886  const Point mapped_point = FEInterface::inverse_map(Dim-1, fe_type,
887  parent_side.get(),
888  support_point);
889 
890  // Compute the parent's side shape function values.
891  for (unsigned int their_side_n=0;
892  their_side_n < n_side_nodes;
893  their_side_n++)
894  {
895  libmesh_assert_less (their_side_n, FEInterface::n_dofs(Dim-1, fe_type, parent_side->type()));
896 
897  const Node * their_node = parent_nodes[their_side_n];
898  libmesh_assert(their_node);
899 
900  const Real their_value = FEInterface::shape(Dim-1,
901  fe_type,
902  parent_side->type(),
903  their_side_n,
904  mapped_point);
905 
906  const Real their_mag = std::abs(their_value);
907 #ifdef DEBUG
908  // Protect for the case u_i ~= u_j,
909  // in which case i better equal j.
910  if (their_mag > 0.999)
911  {
912  libmesh_assert_equal_to (my_node, their_node);
913  libmesh_assert_less (std::abs(their_value - 1.), 0.001);
914  }
915  else
916 #endif
917  // To make nodal constraints useful for constructing
918  // sparsity patterns faster, we need to get EVERY
919  // POSSIBLE constraint coupling identified, even if
920  // there is no coupling in the isoparametric
921  // Lagrange case.
922  if (their_mag < 1.e-5)
923  {
924  // since we may be running this method concurrently
925  // on multiple threads we need to acquire a lock
926  // before modifying the shared constraint_row object.
927  Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
928 
929  // A reference to the constraint row.
930  NodeConstraintRow & constraint_row = constraints[my_node].first;
931 
932  constraint_row.insert(std::make_pair (their_node,
933  0.));
934  }
935  // To get nodal coordinate constraints right, only
936  // add non-zero and non-identity values for Lagrange
937  // basis functions.
938  else // (1.e-5 <= their_mag <= .999)
939  {
940  // since we may be running this method concurrently
941  // on multiple threads we need to acquire a lock
942  // before modifying the shared constraint_row object.
943  Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
944 
945  // A reference to the constraint row.
946  NodeConstraintRow & constraint_row = constraints[my_node].first;
947 
948  constraint_row.insert(std::make_pair (their_node,
949  their_value));
950  }
951  }
952  }
953  }
954 }
955 
956 #endif // LIBMESH_ENABLE_NODE_CONSTRAINTS
957 
958 #endif // #ifdef LIBMESH_ENABLE_AMR
959 
960 
961 
962 #ifdef LIBMESH_ENABLE_PERIODIC
963 
964 #ifdef LIBMESH_ENABLE_NODE_CONSTRAINTS
965 void FEAbstract::compute_periodic_node_constraints (NodeConstraints & constraints,
966  const PeriodicBoundaries & boundaries,
967  const MeshBase & mesh,
968  const PointLocatorBase * point_locator,
969  const Elem * elem)
970 {
971  // Only bother if we truly have periodic boundaries
972  if (boundaries.empty())
973  return;
974 
975  libmesh_assert(elem);
976 
977  // Only constrain active elements with this method
978  if (!elem->active())
979  return;
980 
981  const unsigned int Dim = elem->dim();
982 
983  // We currently always use LAGRANGE mappings for geometry
984  const FEType fe_type(elem->default_order(), LAGRANGE);
985 
986  // Pull objects out of the loop to reduce heap operations
987  std::vector<const Node *> my_nodes, neigh_nodes;
988  std::unique_ptr<const Elem> my_side, neigh_side;
989 
990  // Look at the element faces. Check to see if we need to
991  // build constraints.
992  std::vector<boundary_id_type> bc_ids;
993  for (auto s : elem->side_index_range())
994  {
995  if (elem->neighbor_ptr(s))
996  continue;
997 
998  mesh.get_boundary_info().boundary_ids (elem, s, bc_ids);
999  for (const auto & boundary_id : bc_ids)
1000  {
1001  const PeriodicBoundaryBase * periodic = boundaries.boundary(boundary_id);
1002  if (periodic)
1003  {
1004  libmesh_assert(point_locator);
1005 
1006  // Get pointers to the element's neighbor.
1007  const Elem * neigh = boundaries.neighbor(boundary_id, *point_locator, elem, s);
1008 
1009  // h refinement constraints:
1010  // constrain dofs shared between
1011  // this element and ones as coarse
1012  // as or coarser than this element.
1013  if (neigh->level() <= elem->level())
1014  {
1015  unsigned int s_neigh =
1016  mesh.get_boundary_info().side_with_boundary_id(neigh, periodic->pairedboundary);
1017  libmesh_assert_not_equal_to (s_neigh, libMesh::invalid_uint);
1018 
1019 #ifdef LIBMESH_ENABLE_AMR
1020  libmesh_assert(neigh->active());
1021 #endif // #ifdef LIBMESH_ENABLE_AMR
1022 
1023  elem->build_side_ptr(my_side, s);
1024  neigh->build_side_ptr(neigh_side, s_neigh);
1025 
1026  const unsigned int n_side_nodes = my_side->n_nodes();
1027 
1028  my_nodes.clear();
1029  my_nodes.reserve (n_side_nodes);
1030  neigh_nodes.clear();
1031  neigh_nodes.reserve (n_side_nodes);
1032 
1033  for (unsigned int n=0; n != n_side_nodes; ++n)
1034  my_nodes.push_back(my_side->node_ptr(n));
1035 
1036  for (unsigned int n=0; n != n_side_nodes; ++n)
1037  neigh_nodes.push_back(neigh_side->node_ptr(n));
1038 
1039  // Make sure we're not adding recursive constraints
1040  // due to the redundancy in the way we add periodic
1041  // boundary constraints, or adding constraints to
1042  // nodes that already have AMR constraints
1043  std::vector<bool> skip_constraint(n_side_nodes, false);
1044 
1045  for (unsigned int my_side_n=0;
1046  my_side_n < n_side_nodes;
1047  my_side_n++)
1048  {
1049  libmesh_assert_less (my_side_n, FEInterface::n_dofs(Dim-1, fe_type, my_side->type()));
1050 
1051  const Node * my_node = my_nodes[my_side_n];
1052 
1053  // Figure out where my node lies on their reference element.
1054  const Point neigh_point = periodic->get_corresponding_pos(*my_node);
1055 
1056  const Point mapped_point = FEInterface::inverse_map(Dim-1, fe_type,
1057  neigh_side.get(),
1058  neigh_point);
1059 
1060  // If we've already got a constraint on this
1061  // node, then the periodic constraint is
1062  // redundant
1063  {
1064  Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
1065 
1066  if (constraints.count(my_node))
1067  {
1068  skip_constraint[my_side_n] = true;
1069  continue;
1070  }
1071  }
1072 
1073  // Compute the neighbors's side shape function values.
1074  for (unsigned int their_side_n=0;
1075  their_side_n < n_side_nodes;
1076  their_side_n++)
1077  {
1078  libmesh_assert_less (their_side_n, FEInterface::n_dofs(Dim-1, fe_type, neigh_side->type()));
1079 
1080  const Node * their_node = neigh_nodes[their_side_n];
1081 
1082  // If there's a constraint on an opposing node,
1083  // we need to see if it's constrained by
1084  // *our side* making any periodic constraint
1085  // on us recursive
1086  {
1087  Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
1088 
1089  if (!constraints.count(their_node))
1090  continue;
1091 
1092  const NodeConstraintRow & their_constraint_row =
1093  constraints[their_node].first;
1094 
1095  for (unsigned int orig_side_n=0;
1096  orig_side_n < n_side_nodes;
1097  orig_side_n++)
1098  {
1099  libmesh_assert_less (orig_side_n, FEInterface::n_dofs(Dim-1, fe_type, my_side->type()));
1100 
1101  const Node * orig_node = my_nodes[orig_side_n];
1102 
1103  if (their_constraint_row.count(orig_node))
1104  skip_constraint[orig_side_n] = true;
1105  }
1106  }
1107  }
1108  }
1109  for (unsigned int my_side_n=0;
1110  my_side_n < n_side_nodes;
1111  my_side_n++)
1112  {
1113  libmesh_assert_less (my_side_n, FEInterface::n_dofs(Dim-1, fe_type, my_side->type()));
1114 
1115  if (skip_constraint[my_side_n])
1116  continue;
1117 
1118  const Node * my_node = my_nodes[my_side_n];
1119 
1120  // Figure out where my node lies on their reference element.
1121  const Point neigh_point = periodic->get_corresponding_pos(*my_node);
1122 
1123  // Figure out where my node lies on their reference element.
1124  const Point mapped_point = FEInterface::inverse_map(Dim-1, fe_type,
1125  neigh_side.get(),
1126  neigh_point);
1127 
1128  for (unsigned int their_side_n=0;
1129  their_side_n < n_side_nodes;
1130  their_side_n++)
1131  {
1132  libmesh_assert_less (their_side_n, FEInterface::n_dofs(Dim-1, fe_type, neigh_side->type()));
1133 
1134  const Node * their_node = neigh_nodes[their_side_n];
1135  libmesh_assert(their_node);
1136 
1137  const Real their_value = FEInterface::shape(Dim-1,
1138  fe_type,
1139  neigh_side->type(),
1140  their_side_n,
1141  mapped_point);
1142 
1143  // since we may be running this method concurrently
1144  // on multiple threads we need to acquire a lock
1145  // before modifying the shared constraint_row object.
1146  {
1147  Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
1148 
1149  NodeConstraintRow & constraint_row =
1150  constraints[my_node].first;
1151 
1152  constraint_row.insert(std::make_pair(their_node,
1153  their_value));
1154  }
1155  }
1156  }
1157  }
1158  }
1159  }
1160  }
1161 }
1162 #endif // LIBMESH_ENABLE_NODE_CONSTRAINTS
1163 
1164 #endif // LIBMESH_ENABLE_PERIODIC
1165 
1166 
1167 } // namespace libMesh
libMesh::FEType
Manages the family, order, etc. parameters for a given FE.
Definition: fe_type.h:179
libMesh::FEType::family
FEFamily family
Definition: fe_type.h:204
libMesh::FEAbstract::print_phi
virtual void print_phi(std::ostream &os) const =0
std::abs
double abs(double a)
Definition: libmesh_augment_std_namespace.h:59
libMesh::Elem::parent
const Elem * parent() const
Definition: elem.h:2479
periodic_boundary.h
libMesh::Node
A geometric point in (x,y,z) space associated with a DOF.
Definition: node.h:52
libMesh::FEInterface::n_dofs
static unsigned int n_dofs(const unsigned int dim, const FEType &fe_t, const ElemType t)
Definition: fe_interface.C:454
libMesh::invalid_uint
const unsigned int invalid_uint
Definition: libmesh.h:245
libMesh::FEAbstract::compute_node_constraints
static void compute_node_constraints(NodeConstraints &constraints, const Elem *elem)
Definition: fe_abstract.C:820
libMesh::Elem::side_index_range
IntRange< unsigned short > side_index_range() const
Definition: elem.h:2166
mesh_base.h
libMesh::PeriodicBoundaries
Maps between boundary ids and PeriodicBoundaryBase objects.
Definition: periodic_boundaries.h:51
libMesh::Elem
The base class for all geometric element types.
Definition: elem.h:100
mesh
MeshBase & mesh
Definition: mesh_communication.C:1245
libMesh::PeriodicBoundaries::boundary
PeriodicBoundaryBase * boundary(boundary_id_type id)
Definition: periodic_boundaries.C:40
libMesh::FEAbstract::FEAbstract
FEAbstract(const unsigned int dim, const FEType &fet)
Definition: fe_abstract.C:46
remote_elem.h
fe_interface.h
libMesh::FEAbstract::on_reference_element
static bool on_reference_element(const Point &p, const ElemType t, const Real eps=TOLERANCE)
Definition: fe_abstract.C:581
libMesh::MeshBase
Base class for Mesh.
Definition: mesh_base.h:77
tensor_value.h
dense_matrix.h
libMesh::Elem::build_side_ptr
virtual std::unique_ptr< Elem > build_side_ptr(const unsigned int i, bool proxy=true)=0
libMesh::PeriodicBoundaryBase::get_corresponding_pos
virtual Point get_corresponding_pos(const Point &pt) const =0
libMesh::FEAbstract::print_dphi
virtual void print_dphi(std::ostream &os) const =0
libMesh::Threads::spin_mtx
spin_mutex spin_mtx
Definition: threads.C:29
libMesh::FEAbstract::dim
const unsigned int dim
Definition: fe_abstract.h:531
libMesh::FEAbstract::print_xyz
void print_xyz(std::ostream &os) const
Definition: fe_abstract.C:787
enum_elem_type.h
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::FEAbstract::get_refspace_nodes
static void get_refspace_nodes(const ElemType t, std::vector< Point > &nodes)
Definition: fe_abstract.C:283
libmesh_logging.h
libMesh::FEInterface::inverse_map
static Point inverse_map(const unsigned int dim, const FEType &fe_t, const Elem *elem, const Point &p, const Real tolerance=TOLERANCE, const bool secure=true)
Definition: fe_interface.C:590
periodic_boundaries.h
libMesh::FEAbstract::print_JxW
void print_JxW(std::ostream &os) const
Definition: fe_abstract.C:780
libMesh::FEAbstract::print_info
void print_info(std::ostream &os) const
Definition: fe_abstract.C:793
dof_map.h
libMesh::Elem::neighbor_ptr
const Elem * neighbor_ptr(unsigned int i) const
Definition: elem.h:2050
libMesh::Elem::level
unsigned int level() const
Definition: elem.h:2521
libMesh::Real
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
Definition: libmesh_common.h:131
libMesh::PeriodicBoundaries::neighbor
const Elem * neighbor(boundary_id_type boundary_id, const PointLocatorBase &point_locator, const Elem *e, unsigned int side) const
Definition: periodic_boundaries.C:61
libMesh::Elem::dim
virtual unsigned short dim() const =0
libMesh::Elem::subactive
bool subactive() const
Definition: elem.h:2408
libMesh::FEAbstract::build
static std::unique_ptr< FEAbstract > build(const unsigned int dim, const FEType &type)
Definition: fe_abstract.C:71
libMesh::FEAbstract::~FEAbstract
virtual ~FEAbstract()
Definition: fe_abstract.C:66
boundary_info.h
libMesh::PeriodicBoundaryBase
Base class for all PeriodicBoundary implementations.
Definition: periodic_boundary_base.h:49
libMesh::NodeConstraintRow
std::map< const Node *, Real, std::less< const Node * >, Threads::scalable_allocator< std::pair< const Node *const, Real > > > NodeConstraintRow
Definition: dof_map.h:145
quadrature_gauss.h
dense_vector.h
libMesh::FEAbstract::_fe_map
std::unique_ptr< FEMap > _fe_map
Definition: fe_abstract.h:525
libMesh::Elem::default_order
virtual Order default_order() const =0
libMesh::Elem::active
bool active() const
Definition: elem.h:2390
libMesh::FEMap
Computes finite element mapping function values, gradients, etc.
Definition: fe_map.h:48
libMesh::operator<<
std::ostream & operator<<(std::ostream &os, const FEAbstract &fe)
Definition: fe_abstract.C:809
numeric_vector.h
libMesh::Point
A geometric point in (x,y,z) space.
Definition: point.h:38
libMesh::FEAbstract::compute_periodic_node_constraints
static void compute_periodic_node_constraints(NodeConstraints &constraints, const PeriodicBoundaries &boundaries, const MeshBase &mesh, const PointLocatorBase *point_locator, const Elem *elem)
Definition: fe_abstract.C:965
libMesh::remote_elem
const RemoteElem * remote_elem
Definition: remote_elem.C:57
generated by doxygen