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

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

Hosted By:
SourceForge.net Logo