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