fe_base.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"
22 #include "libmesh/inf_fe.h"
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/tensor_value.h"
38 #include "libmesh/threads.h"
39 
40 // Anonymous namespace, for a helper function for periodic boundary
41 // constraint calculations
42 namespace
43 {
44  using namespace libMesh;
45 
46  // Find the "primary" element around a boundary point:
47  const Elem* primary_boundary_point_neighbor
48  (const Elem* elem,
49  const Point& p,
50  const BoundaryInfo& boundary_info,
51  const std::set<boundary_id_type>& boundary_ids)
52  {
53  // If we don't find a better alternative, the user will have
54  // provided the primary element
55  const Elem *primary = elem;
56 
57  std::set<const Elem*> point_neighbors;
58  elem->find_point_neighbors(p, point_neighbors);
59  for (std::set<const Elem*>::const_iterator point_neighbors_iter =
60  point_neighbors.begin();
61  point_neighbors_iter != point_neighbors.end(); ++point_neighbors_iter)
62  {
63  const Elem* pt_neighbor = *point_neighbors_iter;
64 
65  // If this point neighbor isn't at least
66  // as coarse as the current primary elem, or if it is at
67  // the same level but has a lower id, then
68  // we won't defer to it.
69  if ((pt_neighbor->level() > primary->level()) ||
70  (pt_neighbor->level() == primary->level() &&
71  pt_neighbor->id() < primary->id()))
72  continue;
73 
74  // Otherwise, we will defer to the point neighbor, but only if
75  // one of its sides is on a relevant boundary and that side
76  // contains this vertex
77  bool vertex_on_periodic_side = false;
78  for (unsigned int ns = 0;
79  ns != pt_neighbor->n_sides(); ++ns)
80  {
81  const std::vector<boundary_id_type> bc_ids =
82  boundary_info.boundary_ids (pt_neighbor, ns);
83 
84  bool on_relevant_boundary = false;
85  for (std::set<boundary_id_type>::const_iterator i =
86  boundary_ids.begin(); i != boundary_ids.end(); ++i)
87  if (std::find(bc_ids.begin(), bc_ids.end(), *i)
88  != bc_ids.end())
89  on_relevant_boundary = true;
90 
91  if (!on_relevant_boundary)
92  continue;
93 
94  if (!pt_neighbor->build_side(ns)->contains_point(p))
95  continue;
96 
97  vertex_on_periodic_side = true;
98  break;
99  }
100 
101  if (vertex_on_periodic_side)
102  primary = pt_neighbor;
103  }
104 
105  return primary;
106  }
107 
108  // Find the "primary" element around a boundary edge:
109  const Elem* primary_boundary_edge_neighbor
110  (const Elem* elem,
111  const Point& p1,
112  const Point& p2,
113  const BoundaryInfo& boundary_info,
114  const std::set<boundary_id_type>& boundary_ids)
115  {
116  // If we don't find a better alternative, the user will have
117  // provided the primary element
118  const Elem *primary = elem;
119 
120  std::set<const Elem*> edge_neighbors;
121  elem->find_edge_neighbors(p1, p2, edge_neighbors);
122  for (std::set<const Elem*>::const_iterator edge_neighbors_iter =
123  edge_neighbors.begin();
124  edge_neighbors_iter != edge_neighbors.end(); ++edge_neighbors_iter)
125  {
126  const Elem* e_neighbor = *edge_neighbors_iter;
127 
128  // If this edge neighbor isn't at least
129  // as coarse as the current primary elem, or if it is at
130  // the same level but has a lower id, then
131  // we won't defer to it.
132  if ((e_neighbor->level() > primary->level()) ||
133  (e_neighbor->level() == primary->level() &&
134  e_neighbor->id() < primary->id()))
135  continue;
136 
137  // Otherwise, we will defer to the edge neighbor, but only if
138  // one of its sides is on this periodic boundary and that
139  // side contains this edge
140  bool vertex_on_periodic_side = false;
141  for (unsigned int ns = 0;
142  ns != e_neighbor->n_sides(); ++ns)
143  {
144  const std::vector<boundary_id_type>& bc_ids =
145  boundary_info.boundary_ids (e_neighbor, ns);
146 
147  bool on_relevant_boundary = false;
148  for (std::set<boundary_id_type>::const_iterator i =
149  boundary_ids.begin(); i != boundary_ids.end(); ++i)
150  if (std::find(bc_ids.begin(), bc_ids.end(), *i)
151  != bc_ids.end())
152  on_relevant_boundary = true;
153 
154  if (!on_relevant_boundary)
155  continue;
156 
157  AutoPtr<Elem> periodic_side = e_neighbor->build_side(ns);
158  if (!(periodic_side->contains_point(p1) &&
159  periodic_side->contains_point(p2)))
160  continue;
161 
162  vertex_on_periodic_side = true;
163  break;
164  }
165 
166  if (vertex_on_periodic_side)
167  primary = e_neighbor;
168  }
169 
170  return primary;
171  }
172 
173 }
174 
175 namespace libMesh
176 {
177 
178 
179 
180 // ------------------------------------------------------------
181 // FEBase class members
182 template <>
184 FEGenericBase<Real>::build (const unsigned int dim,
185  const FEType& fet)
186 {
187  // The stupid AutoPtr<FEBase> ap(); return ap;
188  // construct is required to satisfy IBM's xlC
189 
190  switch (dim)
191  {
192  // 0D
193  case 0:
194  {
195  switch (fet.family)
196  {
197  case CLOUGH:
198  {
199  AutoPtr<FEBase> ap(new FE<0,CLOUGH>(fet));
200  return ap;
201  }
202 
203  case HERMITE:
204  {
205  AutoPtr<FEBase> ap(new FE<0,HERMITE>(fet));
206  return ap;
207  }
208 
209  case LAGRANGE:
210  {
211  AutoPtr<FEBase> ap(new FE<0,LAGRANGE>(fet));
212  return ap;
213  }
214 
215  case L2_LAGRANGE:
216  {
217  AutoPtr<FEBase> ap(new FE<0,L2_LAGRANGE>(fet));
218  return ap;
219  }
220 
221  case HIERARCHIC:
222  {
223  AutoPtr<FEBase> ap(new FE<0,HIERARCHIC>(fet));
224  return ap;
225  }
226 
227  case L2_HIERARCHIC:
228  {
230  return ap;
231  }
232 
233  case MONOMIAL:
234  {
235  AutoPtr<FEBase> ap(new FE<0,MONOMIAL>(fet));
236  return ap;
237  }
238 
239 #ifdef LIBMESH_ENABLE_HIGHER_ORDER_SHAPES
240  case SZABAB:
241  {
242  AutoPtr<FEBase> ap(new FE<0,SZABAB>(fet));
243  return ap;
244  }
245 
246  case BERNSTEIN:
247  {
248  AutoPtr<FEBase> ap(new FE<0,BERNSTEIN>(fet));
249  return ap;
250  }
251 #endif
252 
253  case XYZ:
254  {
255  AutoPtr<FEBase> ap(new FEXYZ<0>(fet));
256  return ap;
257  }
258 
259  case SCALAR:
260  {
261  AutoPtr<FEBase> ap(new FEScalar<0>(fet));
262  return ap;
263  }
264 
265  default:
266  libMesh::out << "ERROR: Bad FEType.family= " << fet.family << std::endl;
267  libmesh_error();
268  }
269  }
270  // 1D
271  case 1:
272  {
273  switch (fet.family)
274  {
275  case CLOUGH:
276  {
277  AutoPtr<FEBase> ap(new FE<1,CLOUGH>(fet));
278  return ap;
279  }
280 
281  case HERMITE:
282  {
283  AutoPtr<FEBase> ap(new FE<1,HERMITE>(fet));
284  return ap;
285  }
286 
287  case LAGRANGE:
288  {
289  AutoPtr<FEBase> ap(new FE<1,LAGRANGE>(fet));
290  return ap;
291  }
292 
293  case L2_LAGRANGE:
294  {
295  AutoPtr<FEBase> ap(new FE<1,L2_LAGRANGE>(fet));
296  return ap;
297  }
298 
299  case HIERARCHIC:
300  {
301  AutoPtr<FEBase> ap(new FE<1,HIERARCHIC>(fet));
302  return ap;
303  }
304 
305  case L2_HIERARCHIC:
306  {
308  return ap;
309  }
310 
311  case MONOMIAL:
312  {
313  AutoPtr<FEBase> ap(new FE<1,MONOMIAL>(fet));
314  return ap;
315  }
316 
317 #ifdef LIBMESH_ENABLE_HIGHER_ORDER_SHAPES
318  case SZABAB:
319  {
320  AutoPtr<FEBase> ap(new FE<1,SZABAB>(fet));
321  return ap;
322  }
323 
324  case BERNSTEIN:
325  {
326  AutoPtr<FEBase> ap(new FE<1,BERNSTEIN>(fet));
327  return ap;
328  }
329 #endif
330 
331  case XYZ:
332  {
333  AutoPtr<FEBase> ap(new FEXYZ<1>(fet));
334  return ap;
335  }
336 
337  case SCALAR:
338  {
339  AutoPtr<FEBase> ap(new FEScalar<1>(fet));
340  return ap;
341  }
342 
343  default:
344  libMesh::out << "ERROR: Bad FEType.family= " << fet.family << std::endl;
345  libmesh_error();
346  }
347  }
348 
349 
350  // 2D
351  case 2:
352  {
353  switch (fet.family)
354  {
355  case CLOUGH:
356  {
357  AutoPtr<FEBase> ap(new FE<2,CLOUGH>(fet));
358  return ap;
359  }
360 
361  case HERMITE:
362  {
363  AutoPtr<FEBase> ap(new FE<2,HERMITE>(fet));
364  return ap;
365  }
366 
367  case LAGRANGE:
368  {
369  AutoPtr<FEBase> ap(new FE<2,LAGRANGE>(fet));
370  return ap;
371  }
372 
373  case L2_LAGRANGE:
374  {
375  AutoPtr<FEBase> ap(new FE<2,L2_LAGRANGE>(fet));
376  return ap;
377  }
378 
379  case HIERARCHIC:
380  {
381  AutoPtr<FEBase> ap(new FE<2,HIERARCHIC>(fet));
382  return ap;
383  }
384 
385  case L2_HIERARCHIC:
386  {
388  return ap;
389  }
390 
391  case MONOMIAL:
392  {
393  AutoPtr<FEBase> ap(new FE<2,MONOMIAL>(fet));
394  return ap;
395  }
396 
397 #ifdef LIBMESH_ENABLE_HIGHER_ORDER_SHAPES
398  case SZABAB:
399  {
400  AutoPtr<FEBase> ap(new FE<2,SZABAB>(fet));
401  return ap;
402  }
403 
404  case BERNSTEIN:
405  {
406  AutoPtr<FEBase> ap(new FE<2,BERNSTEIN>(fet));
407  return ap;
408  }
409 #endif
410 
411  case XYZ:
412  {
413  AutoPtr<FEBase> ap(new FEXYZ<2>(fet));
414  return ap;
415  }
416 
417  case SCALAR:
418  {
419  AutoPtr<FEBase> ap(new FEScalar<2>(fet));
420  return ap;
421  }
422  default:
423  libMesh::out << "ERROR: Bad FEType.family= " << fet.family << std::endl;
424  libmesh_error();
425  }
426  }
427 
428 
429  // 3D
430  case 3:
431  {
432  switch (fet.family)
433  {
434  case CLOUGH:
435  {
436  libMesh::out << "ERROR: Clough-Tocher elements currently only support 1D and 2D"
437  << std::endl;
438  libmesh_error();
439  }
440 
441  case HERMITE:
442  {
443  AutoPtr<FEBase> ap(new FE<3,HERMITE>(fet));
444  return ap;
445  }
446 
447  case LAGRANGE:
448  {
449  AutoPtr<FEBase> ap(new FE<3,LAGRANGE>(fet));
450  return ap;
451  }
452 
453  case L2_LAGRANGE:
454  {
455  AutoPtr<FEBase> ap(new FE<3,L2_LAGRANGE>(fet));
456  return ap;
457  }
458 
459  case HIERARCHIC:
460  {
461  AutoPtr<FEBase> ap(new FE<3,HIERARCHIC>(fet));
462  return ap;
463  }
464 
465  case L2_HIERARCHIC:
466  {
468  return ap;
469  }
470 
471  case MONOMIAL:
472  {
473  AutoPtr<FEBase> ap(new FE<3,MONOMIAL>(fet));
474  return ap;
475  }
476 
477 #ifdef LIBMESH_ENABLE_HIGHER_ORDER_SHAPES
478  case SZABAB:
479  {
480  AutoPtr<FEBase> ap(new FE<3,SZABAB>(fet));
481  return ap;
482  }
483 
484  case BERNSTEIN:
485  {
486  AutoPtr<FEBase> ap(new FE<3,BERNSTEIN>(fet));
487  return ap;
488  }
489 #endif
490 
491  case XYZ:
492  {
493  AutoPtr<FEBase> ap(new FEXYZ<3>(fet));
494  return ap;
495  }
496 
497  case SCALAR:
498  {
499  AutoPtr<FEBase> ap(new FEScalar<3>(fet));
500  return ap;
501  }
502 
503  default:
504  libMesh::out << "ERROR: Bad FEType.family= " << fet.family << std::endl;
505  libmesh_error();
506  }
507  }
508 
509  default:
510  libmesh_error();
511  }
512 
513  libmesh_error();
514  AutoPtr<FEBase> ap(NULL);
515  return ap;
516 }
517 
518 
519 
520 template <>
523  const FEType& fet)
524 {
525  // The stupid AutoPtr<FEBase> ap(); return ap;
526  // construct is required to satisfy IBM's xlC
527 
528  switch (dim)
529  {
530  // 0D
531  case 0:
532  {
533  switch (fet.family)
534  {
535  case LAGRANGE_VEC:
536  {
538  return ap;
539  }
540  default:
541  {
542  libMesh::out << "ERROR: Bad FEType.family= " << fet.family << std::endl;
543  libmesh_error();
544  }
545  }
546  }
547  case 1:
548  {
549  switch (fet.family)
550  {
551  case LAGRANGE_VEC:
552  {
554  return ap;
555  }
556  default:
557  {
558  libMesh::out << "ERROR: Bad FEType.family= " << fet.family << std::endl;
559  libmesh_error();
560  }
561  }
562  }
563  case 2:
564  {
565  switch (fet.family)
566  {
567  case LAGRANGE_VEC:
568  {
570  return ap;
571  }
572  case NEDELEC_ONE:
573  {
574  AutoPtr<FEVectorBase> ap( new FENedelecOne<2>(fet) );
575  return ap;
576  }
577  default:
578  {
579  libMesh::out << "ERROR: Bad FEType.family= " << fet.family << std::endl;
580  libmesh_error();
581  }
582  }
583  }
584  case 3:
585  {
586  switch (fet.family)
587  {
588  case LAGRANGE_VEC:
589  {
591  return ap;
592  }
593  case NEDELEC_ONE:
594  {
595  AutoPtr<FEVectorBase> ap( new FENedelecOne<3>(fet) );
596  return ap;
597  }
598  default:
599  {
600  libMesh::out << "ERROR: Bad FEType.family= " << fet.family << std::endl;
601  libmesh_error();
602  }
603  }
604  }
605 
606  default:
607  libmesh_error();
608 
609  } // switch(dim)
610 
611  libmesh_error();
612  AutoPtr<FEVectorBase> ap(NULL);
613  return ap;
614 }
615 
616 
617 
618 
619 
620 
621 
622 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
623 
624 
625 template <>
628  const FEType& fet)
629 {
630  // The stupid AutoPtr<FEBase> ap(); return ap;
631  // construct is required to satisfy IBM's xlC
632 
633  switch (dim)
634  {
635 
636  // 1D
637  case 1:
638  {
639  switch (fet.radial_family)
640  {
641  case INFINITE_MAP:
642  {
643  libMesh::err << "ERROR: Don't build an infinite element " << std::endl
644  << " with FEFamily = " << fet.radial_family << std::endl;
645  libmesh_error();
646  }
647 
648  case JACOBI_20_00:
649  {
650  switch (fet.inf_map)
651  {
652  case CARTESIAN:
653  {
655  return ap;
656  }
657  default:
658  libMesh::err << "ERROR: Don't build an infinite element " << std::endl
659  << " with InfMapType = " << fet.inf_map << std::endl;
660  libmesh_error();
661  }
662  }
663 
664  case JACOBI_30_00:
665  {
666  switch (fet.inf_map)
667  {
668  case CARTESIAN:
669  {
671  return ap;
672  }
673  default:
674  libMesh::err << "ERROR: Don't build an infinite element " << std::endl
675  << " with InfMapType = " << fet.inf_map << std::endl;
676  libmesh_error();
677  }
678  }
679 
680  case LEGENDRE:
681  {
682  switch (fet.inf_map)
683  {
684  case CARTESIAN:
685  {
687  return ap;
688  }
689  default:
690  libMesh::err << "ERROR: Don't build an infinite element " << std::endl
691  << " with InfMapType = " << fet.inf_map << std::endl;
692  libmesh_error();
693  }
694  }
695 
696  case LAGRANGE:
697  {
698  switch (fet.inf_map)
699  {
700  case CARTESIAN:
701  {
703  return ap;
704  }
705  default:
706  libMesh::err << "ERROR: Don't build an infinite element " << std::endl
707  << " with InfMapType = " << fet.inf_map << std::endl;
708  libmesh_error();
709  }
710  }
711 
712 
713 
714  default:
715  libMesh::err << "ERROR: Bad FEType.radial_family= " << fet.radial_family << std::endl;
716  libmesh_error();
717  }
718 
719  }
720 
721 
722 
723 
724  // 2D
725  case 2:
726  {
727  switch (fet.radial_family)
728  {
729  case INFINITE_MAP:
730  {
731  libMesh::err << "ERROR: Don't build an infinite element " << std::endl
732  << " with FEFamily = " << fet.radial_family << std::endl;
733  libmesh_error();
734  }
735 
736  case JACOBI_20_00:
737  {
738  switch (fet.inf_map)
739  {
740  case CARTESIAN:
741  {
743  return ap;
744  }
745  default:
746  libMesh::err << "ERROR: Don't build an infinite element " << std::endl
747  << " with InfMapType = " << fet.inf_map << std::endl;
748  libmesh_error();
749  }
750  }
751 
752  case JACOBI_30_00:
753  {
754  switch (fet.inf_map)
755  {
756  case CARTESIAN:
757  {
759  return ap;
760  }
761  default:
762  libMesh::err << "ERROR: Don't build an infinite element " << std::endl
763  << " with InfMapType = " << fet.inf_map << std::endl;
764  libmesh_error();
765  }
766  }
767 
768  case LEGENDRE:
769  {
770  switch (fet.inf_map)
771  {
772  case CARTESIAN:
773  {
775  return ap;
776  }
777  default:
778  libMesh::err << "ERROR: Don't build an infinite element " << std::endl
779  << " with InfMapType = " << fet.inf_map << std::endl;
780  libmesh_error();
781  }
782  }
783 
784  case LAGRANGE:
785  {
786  switch (fet.inf_map)
787  {
788  case CARTESIAN:
789  {
791  return ap;
792  }
793  default:
794  libMesh::err << "ERROR: Don't build an infinite element " << std::endl
795  << " with InfMapType = " << fet.inf_map << std::endl;
796  libmesh_error();
797  }
798  }
799 
800 
801 
802  default:
803  libMesh::err << "ERROR: Bad FEType.radial_family= " << fet.radial_family << std::endl;
804  libmesh_error();
805  }
806 
807  }
808 
809 
810 
811 
812  // 3D
813  case 3:
814  {
815  switch (fet.radial_family)
816  {
817  case INFINITE_MAP:
818  {
819  libMesh::err << "ERROR: Don't build an infinite element " << std::endl
820  << " with FEFamily = " << fet.radial_family << std::endl;
821  libmesh_error();
822  }
823 
824  case JACOBI_20_00:
825  {
826  switch (fet.inf_map)
827  {
828  case CARTESIAN:
829  {
831  return ap;
832  }
833  default:
834  libMesh::err << "ERROR: Don't build an infinite element " << std::endl
835  << " with InfMapType = " << fet.inf_map << std::endl;
836  libmesh_error();
837  }
838  }
839 
840  case JACOBI_30_00:
841  {
842  switch (fet.inf_map)
843  {
844  case CARTESIAN:
845  {
847  return ap;
848  }
849  default:
850  libMesh::err << "ERROR: Don't build an infinite element " << std::endl
851  << " with InfMapType = " << fet.inf_map << std::endl;
852  libmesh_error();
853  }
854  }
855 
856  case LEGENDRE:
857  {
858  switch (fet.inf_map)
859  {
860  case CARTESIAN:
861  {
863  return ap;
864  }
865  default:
866  libMesh::err << "ERROR: Don't build an infinite element " << std::endl
867  << " with InfMapType = " << fet.inf_map << std::endl;
868  libmesh_error();
869  }
870  }
871 
872  case LAGRANGE:
873  {
874  switch (fet.inf_map)
875  {
876  case CARTESIAN:
877  {
879  return ap;
880  }
881  default:
882  libMesh::err << "ERROR: Don't build an infinite element " << std::endl
883  << " with InfMapType = " << fet.inf_map << std::endl;
884  libmesh_error();
885  }
886  }
887 
888 
889 
890  default:
891  libMesh::err << "ERROR: Bad FEType.radial_family= " << fet.radial_family << std::endl;
892  libmesh_error();
893  }
894  }
895 
896  default:
897  libmesh_error();
898  }
899 
900  libmesh_error();
901  AutoPtr<FEBase> ap(NULL);
902  return ap;
903 }
904 
905 
906 
907 template <>
910  const FEType& )
911 {
912  // No vector types defined... YET.
913  libmesh_error();
914  AutoPtr<FEVectorBase> ap(NULL);
915  return ap;
916 }
917 
918 #endif // ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
919 
920 
921 template <typename OutputType>
923  const std::vector<Point>& qp)
924 {
925  //-------------------------------------------------------------------------
926  // Compute the shape function values (and derivatives)
927  // at the Quadrature points. Note that the actual values
928  // have already been computed via init_shape_functions
929 
930  // Start logging the shape function computation
931  START_LOG("compute_shape_functions()", "FE");
932 
933  calculations_started = true;
934 
935  // If the user forgot to request anything, we'll be safe and
936  // calculate everything:
937 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
938  if (!calculate_phi && !calculate_dphi && !calculate_d2phi && !calculate_curl_phi && !calculate_div_phi)
939  {
940  calculate_phi = calculate_dphi = calculate_d2phi = true;
941  // Only compute curl, div for vector-valued elements
943  {
944  calculate_curl_phi = true;
945  calculate_div_phi = true;
946  }
947  }
948 #else
949  if (!calculate_phi && !calculate_dphi && !calculate_curl_phi && !calculate_div_phi)
950  {
951  calculate_phi = calculate_dphi = true;
952  // Only compute curl for vector-valued elements
954  {
955  calculate_curl_phi = true;
956  calculate_div_phi = true;
957  }
958  }
959 #endif // LIBMESH_ENABLE_SECOND_DERIVATIVES
960 
961 
962  if( calculate_phi )
963  this->_fe_trans->map_phi( this->dim, elem, qp, (*this), this->phi );
964 
965  if( calculate_dphi )
966  this->_fe_trans->map_dphi( this->dim, elem, qp, (*this), this->dphi,
967  this->dphidx, this->dphidy, this->dphidz );
968 
969 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
970  if( calculate_d2phi )
971  this->_fe_trans->map_d2phi( this->dim, elem, qp, (*this), this->d2phi,
972  this->d2phidx2, this->d2phidxdy, this->d2phidxdz,
973  this->d2phidy2, this->d2phidydz, this->d2phidz2 );
974 #endif //LIBMESH_ENABLE_SECOND_DERIVATIVES
975 
976  // Only compute curl for vector-valued elements
977  if( calculate_curl_phi && TypesEqual<OutputType,RealGradient>::value )
978  this->_fe_trans->map_curl( this->dim, elem, qp, (*this), this->curl_phi );
979 
980  // Only compute div for vector-valued elements
981  if( calculate_div_phi && TypesEqual<OutputType,RealGradient>::value )
982  this->_fe_trans->map_div( this->dim, elem, qp, (*this), this->div_phi );
983 
984  // Stop logging the shape function computation
985  STOP_LOG("compute_shape_functions()", "FE");
986 }
987 
988 
989 template <typename OutputType>
990 void FEGenericBase<OutputType>::print_phi(std::ostream& os) const
991 {
992  for (unsigned int i=0; i<phi.size(); ++i)
993  for (unsigned int j=0; j<phi[i].size(); ++j)
994  os << " phi[" << i << "][" << j << "]=" << phi[i][j] << std::endl;
995 }
996 
997 
998 
999 
1000 template <typename OutputType>
1001 void FEGenericBase<OutputType>::print_dphi(std::ostream& os) const
1002 {
1003  for (unsigned int i=0; i<dphi.size(); ++i)
1004  for (unsigned int j=0; j<dphi[i].size(); ++j)
1005  os << " dphi[" << i << "][" << j << "]=" << dphi[i][j];
1006 }
1007 
1008 
1009 
1010 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
1011 
1012 
1013 template <typename OutputType>
1014 void FEGenericBase<OutputType>::print_d2phi(std::ostream& os) const
1015 {
1016  for (unsigned int i=0; i<dphi.size(); ++i)
1017  for (unsigned int j=0; j<dphi[i].size(); ++j)
1018  os << " d2phi[" << i << "][" << j << "]=" << d2phi[i][j];
1019 }
1020 
1021 #endif
1022 
1023 
1024 
1025 #ifdef LIBMESH_ENABLE_AMR
1026 
1027 template <typename OutputType>
1028 void
1030  const DofMap &dof_map,
1031  const Elem *elem,
1032  DenseVector<Number> &Ue,
1033  const unsigned int var,
1034  const bool use_old_dof_indices)
1035 {
1036  // Side/edge local DOF indices
1037  std::vector<unsigned int> new_side_dofs, old_side_dofs;
1038 
1039  // FIXME: what about 2D shells in 3D space?
1040  unsigned int dim = elem->dim();
1041 
1042  // We use local FE objects for now
1043  // FIXME: we should use more, external objects instead for efficiency
1044  const FEType& base_fe_type = dof_map.variable_type(var);
1046  (FEGenericBase<OutputShape>::build(dim, base_fe_type));
1048  (FEGenericBase<OutputShape>::build(dim, base_fe_type));
1049 
1050  AutoPtr<QBase> qrule (base_fe_type.default_quadrature_rule(dim));
1051  AutoPtr<QBase> qedgerule (base_fe_type.default_quadrature_rule(1));
1052  AutoPtr<QBase> qsiderule (base_fe_type.default_quadrature_rule(dim-1));
1053  std::vector<Point> coarse_qpoints;
1054 
1055  // The values of the shape functions at the quadrature
1056  // points
1057  const std::vector<std::vector<OutputShape> >& phi_values =
1058  fe->get_phi();
1059  const std::vector<std::vector<OutputShape> >& phi_coarse =
1060  fe_coarse->get_phi();
1061 
1062  // The gradients of the shape functions at the quadrature
1063  // points on the child element.
1064  const std::vector<std::vector<OutputGradient> > *dphi_values =
1065  NULL;
1066  const std::vector<std::vector<OutputGradient> > *dphi_coarse =
1067  NULL;
1068 
1069  const FEContinuity cont = fe->get_continuity();
1070 
1071  if (cont == C_ONE)
1072  {
1073  const std::vector<std::vector<OutputGradient> >&
1074  ref_dphi_values = fe->get_dphi();
1075  dphi_values = &ref_dphi_values;
1076  const std::vector<std::vector<OutputGradient> >&
1077  ref_dphi_coarse = fe_coarse->get_dphi();
1078  dphi_coarse = &ref_dphi_coarse;
1079  }
1080 
1081  // The Jacobian * quadrature weight at the quadrature points
1082  const std::vector<Real>& JxW =
1083  fe->get_JxW();
1084 
1085  // The XYZ locations of the quadrature points on the
1086  // child element
1087  const std::vector<Point>& xyz_values =
1088  fe->get_xyz();
1089 
1090 
1091 
1092  FEType fe_type = base_fe_type, temp_fe_type;
1093  const ElemType elem_type = elem->type();
1094  fe_type.order = static_cast<Order>(fe_type.order +
1095  elem->max_descendant_p_level());
1096 
1097  // Number of nodes on parent element
1098  const unsigned int n_nodes = elem->n_nodes();
1099 
1100  // Number of dofs on parent element
1101  const unsigned int new_n_dofs =
1102  FEInterface::n_dofs(dim, fe_type, elem_type);
1103 
1104  // Fixed vs. free DoFs on edge/face projections
1105  std::vector<char> dof_is_fixed(new_n_dofs, false); // bools
1106  std::vector<int> free_dof(new_n_dofs, 0);
1107 
1108  DenseMatrix<Real> Ke;
1110  Ue.resize(new_n_dofs); Ue.zero();
1111 
1112 
1113  // When coarsening, in general, we need a series of
1114  // projections to ensure a unique and continuous
1115  // solution. We start by interpolating nodes, then
1116  // hold those fixed and project edges, then
1117  // hold those fixed and project faces, then
1118  // hold those fixed and project interiors
1119 
1120  // Copy node values first
1121  {
1122  std::vector<dof_id_type> node_dof_indices;
1123  if (use_old_dof_indices)
1124  dof_map.old_dof_indices (elem, node_dof_indices, var);
1125  else
1126  dof_map.dof_indices (elem, node_dof_indices, var);
1127 
1128  unsigned int current_dof = 0;
1129  for (unsigned int n=0; n!= n_nodes; ++n)
1130  {
1131  // FIXME: this should go through the DofMap,
1132  // not duplicate dof_indices code badly!
1133  const unsigned int my_nc =
1134  FEInterface::n_dofs_at_node (dim, fe_type,
1135  elem_type, n);
1136  if (!elem->is_vertex(n))
1137  {
1138  current_dof += my_nc;
1139  continue;
1140  }
1141 
1142  temp_fe_type = base_fe_type;
1143  // We're assuming here that child n shares vertex n,
1144  // which is wrong on non-simplices right now
1145  // ... but this code isn't necessary except on elements
1146  // where p refinement creates more vertex dofs; we have
1147  // no such elements yet.
1148 /*
1149  if (elem->child(n)->p_level() < elem->p_level())
1150  {
1151  temp_fe_type.order =
1152  static_cast<Order>(temp_fe_type.order +
1153  elem->child(n)->p_level());
1154  }
1155 */
1156  const unsigned int nc =
1157  FEInterface::n_dofs_at_node (dim, temp_fe_type,
1158  elem_type, n);
1159  for (unsigned int i=0; i!= nc; ++i)
1160  {
1161  Ue(current_dof) =
1162  old_vector(node_dof_indices[current_dof]);
1163  dof_is_fixed[current_dof] = true;
1164  current_dof++;
1165  }
1166  }
1167  }
1168 
1169  // In 3D, project any edge values next
1170  if (dim > 2 && cont != DISCONTINUOUS)
1171  for (unsigned int e=0; e != elem->n_edges(); ++e)
1172  {
1173  FEInterface::dofs_on_edge(elem, dim, fe_type,
1174  e, new_side_dofs);
1175 
1176  // Some edge dofs are on nodes and already
1177  // fixed, others are free to calculate
1178  unsigned int free_dofs = 0;
1179  for (unsigned int i=0; i !=
1180  new_side_dofs.size(); ++i)
1181  if (!dof_is_fixed[new_side_dofs[i]])
1182  free_dof[free_dofs++] = i;
1183  Ke.resize (free_dofs, free_dofs); Ke.zero();
1184  Fe.resize (free_dofs); Fe.zero();
1185  // The new edge coefficients
1186  DenseVector<Number> Uedge(free_dofs);
1187 
1188  // Add projection terms from each child sharing
1189  // this edge
1190  for (unsigned int c=0; c != elem->n_children();
1191  ++c)
1192  {
1193  if (!elem->is_child_on_edge(c,e))
1194  continue;
1195  Elem *child = elem->child(c);
1196 
1197  std::vector<dof_id_type> child_dof_indices;
1198  if (use_old_dof_indices)
1199  dof_map.old_dof_indices (child,
1200  child_dof_indices, var);
1201  else
1202  dof_map.dof_indices (child,
1203  child_dof_indices, var);
1204  const unsigned int child_n_dofs =
1205  libmesh_cast_int<unsigned int>
1206  (child_dof_indices.size());
1207 
1208  temp_fe_type = base_fe_type;
1209  temp_fe_type.order =
1210  static_cast<Order>(temp_fe_type.order +
1211  child->p_level());
1212 
1213  FEInterface::dofs_on_edge(child, dim,
1214  temp_fe_type, e, old_side_dofs);
1215 
1216  // Initialize both child and parent FE data
1217  // on the child's edge
1218  fe->attach_quadrature_rule (qedgerule.get());
1219  fe->edge_reinit (child, e);
1220  const unsigned int n_qp = qedgerule->n_points();
1221 
1222  FEInterface::inverse_map (dim, fe_type, elem,
1223  xyz_values, coarse_qpoints);
1224 
1225  fe_coarse->reinit(elem, &coarse_qpoints);
1226 
1227  // Loop over the quadrature points
1228  for (unsigned int qp=0; qp<n_qp; qp++)
1229  {
1230  // solution value at the quadrature point
1231  OutputNumber fineval = libMesh::zero;
1232  // solution grad at the quadrature point
1233  OutputNumberGradient finegrad;
1234 
1235  // Sum the solution values * the DOF
1236  // values at the quadrature point to
1237  // get the solution value and gradient.
1238  for (unsigned int i=0; i<child_n_dofs;
1239  i++)
1240  {
1241  fineval +=
1242  (old_vector(child_dof_indices[i])*
1243  phi_values[i][qp]);
1244  if (cont == C_ONE)
1245  finegrad += (*dphi_values)[i][qp] *
1246  old_vector(child_dof_indices[i]);
1247  }
1248 
1249  // Form edge projection matrix
1250  for (unsigned int sidei=0, freei=0;
1251  sidei != new_side_dofs.size();
1252  ++sidei)
1253  {
1254  unsigned int i = new_side_dofs[sidei];
1255  // fixed DoFs aren't test functions
1256  if (dof_is_fixed[i])
1257  continue;
1258  for (unsigned int sidej=0, freej=0;
1259  sidej != new_side_dofs.size();
1260  ++sidej)
1261  {
1262  unsigned int j =
1263  new_side_dofs[sidej];
1264  if (dof_is_fixed[j])
1265  Fe(freei) -=
1266  TensorTools::inner_product(phi_coarse[i][qp],
1267  phi_coarse[j][qp]) *
1268  JxW[qp] * Ue(j);
1269  else
1270  Ke(freei,freej) +=
1271  TensorTools::inner_product(phi_coarse[i][qp],
1272  phi_coarse[j][qp]) *
1273  JxW[qp];
1274  if (cont == C_ONE)
1275  {
1276  if (dof_is_fixed[j])
1277  Fe(freei) -=
1278  TensorTools::inner_product((*dphi_coarse)[i][qp],
1279  (*dphi_coarse)[j][qp]) *
1280  JxW[qp] * Ue(j);
1281  else
1282  Ke(freei,freej) +=
1283  TensorTools::inner_product((*dphi_coarse)[i][qp],
1284  (*dphi_coarse)[j][qp]) *
1285  JxW[qp];
1286  }
1287  if (!dof_is_fixed[j])
1288  freej++;
1289  }
1290  Fe(freei) += TensorTools::inner_product(phi_coarse[i][qp],
1291  fineval) * JxW[qp];
1292  if (cont == C_ONE)
1293  Fe(freei) +=
1294  TensorTools::inner_product(finegrad, (*dphi_coarse)[i][qp]) * JxW[qp];
1295  freei++;
1296  }
1297  }
1298  }
1299  Ke.cholesky_solve(Fe, Uedge);
1300 
1301  // Transfer new edge solutions to element
1302  for (unsigned int i=0; i != free_dofs; ++i)
1303  {
1304  Number &ui = Ue(new_side_dofs[free_dof[i]]);
1306  std::abs(ui - Uedge(i)) < TOLERANCE);
1307  ui = Uedge(i);
1308  dof_is_fixed[new_side_dofs[free_dof[i]]] =
1309  true;
1310  }
1311  }
1312 
1313  // Project any side values (edges in 2D, faces in 3D)
1314  if (dim > 1 && cont != DISCONTINUOUS)
1315  for (unsigned int s=0; s != elem->n_sides(); ++s)
1316  {
1317  FEInterface::dofs_on_side(elem, dim, fe_type,
1318  s, new_side_dofs);
1319 
1320  // Some side dofs are on nodes/edges and already
1321  // fixed, others are free to calculate
1322  unsigned int free_dofs = 0;
1323  for (unsigned int i=0; i !=
1324  new_side_dofs.size(); ++i)
1325  if (!dof_is_fixed[new_side_dofs[i]])
1326  free_dof[free_dofs++] = i;
1327  Ke.resize (free_dofs, free_dofs); Ke.zero();
1328  Fe.resize (free_dofs); Fe.zero();
1329  // The new side coefficients
1330  DenseVector<Number> Uside(free_dofs);
1331 
1332  // Add projection terms from each child sharing
1333  // this side
1334  for (unsigned int c=0; c != elem->n_children();
1335  ++c)
1336  {
1337  if (!elem->is_child_on_side(c,s))
1338  continue;
1339  Elem *child = elem->child(c);
1340 
1341  std::vector<dof_id_type> child_dof_indices;
1342  if (use_old_dof_indices)
1343  dof_map.old_dof_indices (child,
1344  child_dof_indices, var);
1345  else
1346  dof_map.dof_indices (child,
1347  child_dof_indices, var);
1348  const unsigned int child_n_dofs =
1349  libmesh_cast_int<unsigned int>
1350  (child_dof_indices.size());
1351 
1352  temp_fe_type = base_fe_type;
1353  temp_fe_type.order =
1354  static_cast<Order>(temp_fe_type.order +
1355  child->p_level());
1356 
1357  FEInterface::dofs_on_side(child, dim,
1358  temp_fe_type, s, old_side_dofs);
1359 
1360  // Initialize both child and parent FE data
1361  // on the child's side
1362  fe->attach_quadrature_rule (qsiderule.get());
1363  fe->reinit (child, s);
1364  const unsigned int n_qp = qsiderule->n_points();
1365 
1366  FEInterface::inverse_map (dim, fe_type, elem,
1367  xyz_values, coarse_qpoints);
1368 
1369  fe_coarse->reinit(elem, &coarse_qpoints);
1370 
1371  // Loop over the quadrature points
1372  for (unsigned int qp=0; qp<n_qp; qp++)
1373  {
1374  // solution value at the quadrature point
1375  OutputNumber fineval = libMesh::zero;
1376  // solution grad at the quadrature point
1377  OutputNumberGradient finegrad;
1378 
1379  // Sum the solution values * the DOF
1380  // values at the quadrature point to
1381  // get the solution value and gradient.
1382  for (unsigned int i=0; i<child_n_dofs;
1383  i++)
1384  {
1385  fineval +=
1386  old_vector(child_dof_indices[i]) *
1387  phi_values[i][qp];
1388  if (cont == C_ONE)
1389  finegrad += (*dphi_values)[i][qp] *
1390  old_vector(child_dof_indices[i]);
1391  }
1392 
1393  // Form side projection matrix
1394  for (unsigned int sidei=0, freei=0;
1395  sidei != new_side_dofs.size();
1396  ++sidei)
1397  {
1398  unsigned int i = new_side_dofs[sidei];
1399  // fixed DoFs aren't test functions
1400  if (dof_is_fixed[i])
1401  continue;
1402  for (unsigned int sidej=0, freej=0;
1403  sidej != new_side_dofs.size();
1404  ++sidej)
1405  {
1406  unsigned int j =
1407  new_side_dofs[sidej];
1408  if (dof_is_fixed[j])
1409  Fe(freei) -=
1410  TensorTools::inner_product(phi_coarse[i][qp],
1411  phi_coarse[j][qp]) *
1412  JxW[qp] * Ue(j);
1413  else
1414  Ke(freei,freej) +=
1415  TensorTools::inner_product(phi_coarse[i][qp],
1416  phi_coarse[j][qp]) *
1417  JxW[qp];
1418  if (cont == C_ONE)
1419  {
1420  if (dof_is_fixed[j])
1421  Fe(freei) -=
1422  TensorTools::inner_product((*dphi_coarse)[i][qp],
1423  (*dphi_coarse)[j][qp]) *
1424  JxW[qp] * Ue(j);
1425  else
1426  Ke(freei,freej) +=
1427  TensorTools::inner_product((*dphi_coarse)[i][qp],
1428  (*dphi_coarse)[j][qp]) *
1429  JxW[qp];
1430  }
1431  if (!dof_is_fixed[j])
1432  freej++;
1433  }
1434  Fe(freei) += TensorTools::inner_product(fineval, phi_coarse[i][qp]) * JxW[qp];
1435  if (cont == C_ONE)
1436  Fe(freei) +=
1437  TensorTools::inner_product(finegrad, (*dphi_coarse)[i][qp]) * JxW[qp];
1438  freei++;
1439  }
1440  }
1441  }
1442  Ke.cholesky_solve(Fe, Uside);
1443 
1444  // Transfer new side solutions to element
1445  for (unsigned int i=0; i != free_dofs; ++i)
1446  {
1447  Number &ui = Ue(new_side_dofs[free_dof[i]]);
1449  std::abs(ui - Uside(i)) < TOLERANCE);
1450  ui = Uside(i);
1451  dof_is_fixed[new_side_dofs[free_dof[i]]] =
1452  true;
1453  }
1454  }
1455 
1456  // Project the interior values, finally
1457 
1458  // Some interior dofs are on nodes/edges/sides and
1459  // already fixed, others are free to calculate
1460  unsigned int free_dofs = 0;
1461  for (unsigned int i=0; i != new_n_dofs; ++i)
1462  if (!dof_is_fixed[i])
1463  free_dof[free_dofs++] = i;
1464  Ke.resize (free_dofs, free_dofs); Ke.zero();
1465  Fe.resize (free_dofs); Fe.zero();
1466  // The new interior coefficients
1467  DenseVector<Number> Uint(free_dofs);
1468 
1469  // Add projection terms from each child
1470  for (unsigned int c=0; c != elem->n_children(); ++c)
1471  {
1472  Elem *child = elem->child(c);
1473 
1474  std::vector<dof_id_type> child_dof_indices;
1475  if (use_old_dof_indices)
1476  dof_map.old_dof_indices (child,
1477  child_dof_indices, var);
1478  else
1479  dof_map.dof_indices (child,
1480  child_dof_indices, var);
1481  const unsigned int child_n_dofs =
1482  libmesh_cast_int<unsigned int>
1483  (child_dof_indices.size());
1484 
1485  // Initialize both child and parent FE data
1486  // on the child's quadrature points
1487  fe->attach_quadrature_rule (qrule.get());
1488  fe->reinit (child);
1489  const unsigned int n_qp = qrule->n_points();
1490 
1491  FEInterface::inverse_map (dim, fe_type, elem,
1492  xyz_values, coarse_qpoints);
1493 
1494  fe_coarse->reinit(elem, &coarse_qpoints);
1495 
1496  // Loop over the quadrature points
1497  for (unsigned int qp=0; qp<n_qp; qp++)
1498  {
1499  // solution value at the quadrature point
1500  OutputNumber fineval = libMesh::zero;
1501  // solution grad at the quadrature point
1502  OutputNumberGradient finegrad;
1503 
1504  // Sum the solution values * the DOF
1505  // values at the quadrature point to
1506  // get the solution value and gradient.
1507  for (unsigned int i=0; i<child_n_dofs; i++)
1508  {
1509  fineval +=
1510  (old_vector(child_dof_indices[i]) *
1511  phi_values[i][qp]);
1512  if (cont == C_ONE)
1513  finegrad += (*dphi_values)[i][qp] *
1514  old_vector(child_dof_indices[i]);
1515  }
1516 
1517  // Form interior projection matrix
1518  for (unsigned int i=0, freei=0;
1519  i != new_n_dofs; ++i)
1520  {
1521  // fixed DoFs aren't test functions
1522  if (dof_is_fixed[i])
1523  continue;
1524  for (unsigned int j=0, freej=0; j !=
1525  new_n_dofs; ++j)
1526  {
1527  if (dof_is_fixed[j])
1528  Fe(freei) -=
1529  TensorTools::inner_product(phi_coarse[i][qp],
1530  phi_coarse[j][qp]) *
1531  JxW[qp] * Ue(j);
1532  else
1533  Ke(freei,freej) +=
1534  TensorTools::inner_product(phi_coarse[i][qp],
1535  phi_coarse[j][qp]) *
1536  JxW[qp];
1537  if (cont == C_ONE)
1538  {
1539  if (dof_is_fixed[j])
1540  Fe(freei) -=
1541  TensorTools::inner_product((*dphi_coarse)[i][qp],
1542  (*dphi_coarse)[j][qp]) *
1543  JxW[qp] * Ue(j);
1544  else
1545  Ke(freei,freej) +=
1546  TensorTools::inner_product((*dphi_coarse)[i][qp],
1547  (*dphi_coarse)[j][qp]) *
1548  JxW[qp];
1549  }
1550  if (!dof_is_fixed[j])
1551  freej++;
1552  }
1553  Fe(freei) += TensorTools::inner_product(phi_coarse[i][qp], fineval) *
1554  JxW[qp];
1555  if (cont == C_ONE)
1556  Fe(freei) += TensorTools::inner_product(finegrad, (*dphi_coarse)[i][qp]) * JxW[qp];
1557  freei++;
1558  }
1559  }
1560  }
1561  Ke.cholesky_solve(Fe, Uint);
1562 
1563  // Transfer new interior solutions to element
1564  for (unsigned int i=0; i != free_dofs; ++i)
1565  {
1566  Number &ui = Ue(free_dof[i]);
1568  std::abs(ui - Uint(i)) < TOLERANCE);
1569  ui = Uint(i);
1570  dof_is_fixed[free_dof[i]] = true;
1571  }
1572 
1573  // Make sure every DoF got reached!
1574  for (unsigned int i=0; i != new_n_dofs; ++i)
1575  libmesh_assert(dof_is_fixed[i]);
1576 }
1577 
1578 
1579 
1580 template <typename OutputType>
1581 void
1583  DofMap &dof_map,
1584  const unsigned int variable_number,
1585  const Elem* elem)
1586 {
1587  libmesh_assert(elem);
1588 
1589  const unsigned int Dim = elem->dim();
1590 
1591  // Only constrain elements in 2,3D.
1592  if (Dim == 1)
1593  return;
1594 
1595  // Only constrain active elements with this method
1596  if (!elem->active())
1597  return;
1598 
1599  const FEType& base_fe_type = dof_map.variable_type(variable_number);
1600 
1601  // Construct FE objects for this element and its neighbors.
1603  (FEGenericBase<OutputShape>::build(Dim, base_fe_type));
1604  const FEContinuity cont = my_fe->get_continuity();
1605 
1606  // We don't need to constrain discontinuous elements
1607  if (cont == DISCONTINUOUS)
1608  return;
1609  libmesh_assert (cont == C_ZERO || cont == C_ONE);
1610 
1612  (FEGenericBase<OutputShape>::build(Dim, base_fe_type));
1613 
1614  QGauss my_qface(Dim-1, base_fe_type.default_quadrature_order());
1615  my_fe->attach_quadrature_rule (&my_qface);
1616  std::vector<Point> neigh_qface;
1617 
1618  const std::vector<Real>& JxW = my_fe->get_JxW();
1619  const std::vector<Point>& q_point = my_fe->get_xyz();
1620  const std::vector<std::vector<OutputShape> >& phi = my_fe->get_phi();
1621  const std::vector<std::vector<OutputShape> >& neigh_phi =
1622  neigh_fe->get_phi();
1623  const std::vector<Point> *face_normals = NULL;
1624  const std::vector<std::vector<OutputGradient> > *dphi = NULL;
1625  const std::vector<std::vector<OutputGradient> > *neigh_dphi = NULL;
1626 
1627  std::vector<dof_id_type> my_dof_indices, neigh_dof_indices;
1628  std::vector<unsigned int> my_side_dofs, neigh_side_dofs;
1629 
1630  if (cont != C_ZERO)
1631  {
1632  const std::vector<Point>& ref_face_normals =
1633  my_fe->get_normals();
1634  face_normals = &ref_face_normals;
1635  const std::vector<std::vector<OutputGradient> >& ref_dphi =
1636  my_fe->get_dphi();
1637  dphi = &ref_dphi;
1638  const std::vector<std::vector<OutputGradient> >& ref_neigh_dphi =
1639  neigh_fe->get_dphi();
1640  neigh_dphi = &ref_neigh_dphi;
1641  }
1642 
1643  DenseMatrix<Real> Ke;
1644  DenseVector<Real> Fe;
1645  std::vector<DenseVector<Real> > Ue;
1646 
1647  // Look at the element faces. Check to see if we need to
1648  // build constraints.
1649  for (unsigned int s=0; s<elem->n_sides(); s++)
1650  if (elem->neighbor(s) != NULL)
1651  {
1652  // Get pointers to the element's neighbor.
1653  const Elem* neigh = elem->neighbor(s);
1654 
1655  // h refinement constraints:
1656  // constrain dofs shared between
1657  // this element and ones coarser
1658  // than this element.
1659  if (neigh->level() < elem->level())
1660  {
1661  unsigned int s_neigh = neigh->which_neighbor_am_i(elem);
1662  libmesh_assert_less (s_neigh, neigh->n_neighbors());
1663 
1664  // Find the minimum p level; we build the h constraint
1665  // matrix with this and then constrain away all higher p
1666  // DoFs.
1667  libmesh_assert(neigh->active());
1668  const unsigned int min_p_level =
1669  std::min(elem->p_level(), neigh->p_level());
1670 
1671  // we may need to make the FE objects reinit with the
1672  // minimum shared p_level
1673  // FIXME - I hate using const_cast<> and avoiding
1674  // accessor functions; there's got to be a
1675  // better way to do this!
1676  const unsigned int old_elem_level = elem->p_level();
1677  if (old_elem_level != min_p_level)
1678  (const_cast<Elem *>(elem))->hack_p_level(min_p_level);
1679  const unsigned int old_neigh_level = neigh->p_level();
1680  if (old_neigh_level != min_p_level)
1681  (const_cast<Elem *>(neigh))->hack_p_level(min_p_level);
1682 
1683  my_fe->reinit(elem, s);
1684 
1685  // This function gets called element-by-element, so there
1686  // will be a lot of memory allocation going on. We can
1687  // at least minimize this for the case of the dof indices
1688  // by efficiently preallocating the requisite storage.
1689  // n_nodes is not necessarily n_dofs, but it is better
1690  // than nothing!
1691  my_dof_indices.reserve (elem->n_nodes());
1692  neigh_dof_indices.reserve (neigh->n_nodes());
1693 
1694  dof_map.dof_indices (elem, my_dof_indices,
1695  variable_number);
1696  dof_map.dof_indices (neigh, neigh_dof_indices,
1697  variable_number);
1698 
1699  const unsigned int n_qp = my_qface.n_points();
1700 
1701  FEInterface::inverse_map (Dim, base_fe_type, neigh,
1702  q_point, neigh_qface);
1703 
1704  neigh_fe->reinit(neigh, &neigh_qface);
1705 
1706  // We're only concerned with DOFs whose values (and/or first
1707  // derivatives for C1 elements) are supported on side nodes
1708  FEInterface::dofs_on_side(elem, Dim, base_fe_type, s, my_side_dofs);
1709  FEInterface::dofs_on_side(neigh, Dim, base_fe_type, s_neigh, neigh_side_dofs);
1710 
1711  // We're done with functions that examine Elem::p_level(),
1712  // so let's unhack those levels
1713  if (elem->p_level() != old_elem_level)
1714  (const_cast<Elem *>(elem))->hack_p_level(old_elem_level);
1715  if (neigh->p_level() != old_neigh_level)
1716  (const_cast<Elem *>(neigh))->hack_p_level(old_neigh_level);
1717 
1718  const unsigned int n_side_dofs =
1719  libmesh_cast_int<unsigned int>(my_side_dofs.size());
1720  libmesh_assert_equal_to (n_side_dofs, neigh_side_dofs.size());
1721 
1722  Ke.resize (n_side_dofs, n_side_dofs);
1723  Ue.resize(n_side_dofs);
1724 
1725  // Form the projection matrix, (inner product of fine basis
1726  // functions against fine test functions)
1727  for (unsigned int is = 0; is != n_side_dofs; ++is)
1728  {
1729  const unsigned int i = my_side_dofs[is];
1730  for (unsigned int js = 0; js != n_side_dofs; ++js)
1731  {
1732  const unsigned int j = my_side_dofs[js];
1733  for (unsigned int qp = 0; qp != n_qp; ++qp)
1734  {
1735  Ke(is,js) += JxW[qp] * TensorTools::inner_product(phi[i][qp], phi[j][qp]);
1736  if (cont != C_ZERO)
1737  Ke(is,js) += JxW[qp] *
1738  TensorTools::inner_product((*dphi)[i][qp] *
1739  (*face_normals)[qp],
1740  (*dphi)[j][qp] *
1741  (*face_normals)[qp]);
1742  }
1743  }
1744  }
1745 
1746  // Form the right hand sides, (inner product of coarse basis
1747  // functions against fine test functions)
1748  for (unsigned int is = 0; is != n_side_dofs; ++is)
1749  {
1750  const unsigned int i = neigh_side_dofs[is];
1751  Fe.resize (n_side_dofs);
1752  for (unsigned int js = 0; js != n_side_dofs; ++js)
1753  {
1754  const unsigned int j = my_side_dofs[js];
1755  for (unsigned int qp = 0; qp != n_qp; ++qp)
1756  {
1757  Fe(js) += JxW[qp] *
1758  TensorTools::inner_product(neigh_phi[i][qp],
1759  phi[j][qp]);
1760  if (cont != C_ZERO)
1761  Fe(js) += JxW[qp] *
1762  TensorTools::inner_product((*neigh_dphi)[i][qp] *
1763  (*face_normals)[qp],
1764  (*dphi)[j][qp] *
1765  (*face_normals)[qp]);
1766  }
1767  }
1768  Ke.cholesky_solve(Fe, Ue[is]);
1769  }
1770 
1771  for (unsigned int js = 0; js != n_side_dofs; ++js)
1772  {
1773  const unsigned int j = my_side_dofs[js];
1774  const dof_id_type my_dof_g = my_dof_indices[j];
1775  libmesh_assert_not_equal_to (my_dof_g, DofObject::invalid_id);
1776 
1777  // Hunt for "constraining against myself" cases before
1778  // we bother creating a constraint row
1779  bool self_constraint = false;
1780  for (unsigned int is = 0; is != n_side_dofs; ++is)
1781  {
1782  const unsigned int i = neigh_side_dofs[is];
1783  const dof_id_type their_dof_g = neigh_dof_indices[i];
1784  libmesh_assert_not_equal_to (their_dof_g, DofObject::invalid_id);
1785 
1786  if (their_dof_g == my_dof_g)
1787  {
1788 #ifndef NDEBUG
1789  const Real their_dof_value = Ue[is](js);
1790  libmesh_assert_less (std::abs(their_dof_value-1.),
1791  10*TOLERANCE);
1792 
1793  for (unsigned int k = 0; k != n_side_dofs; ++k)
1794  libmesh_assert(k == is ||
1795  std::abs(Ue[k](js)) <
1796  10*TOLERANCE);
1797 #endif
1798 
1799  self_constraint = true;
1800  break;
1801  }
1802  }
1803 
1804  if (self_constraint)
1805  continue;
1806 
1807  DofConstraintRow* constraint_row;
1808 
1809  // we may be running constraint methods concurretly
1810  // on multiple threads, so we need a lock to
1811  // ensure that this constraint is "ours"
1812  {
1813  Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
1814 
1815  if (dof_map.is_constrained_dof(my_dof_g))
1816  continue;
1817 
1818  constraint_row = &(constraints[my_dof_g]);
1819  libmesh_assert(constraint_row->empty());
1820  }
1821 
1822  for (unsigned int is = 0; is != n_side_dofs; ++is)
1823  {
1824  const unsigned int i = neigh_side_dofs[is];
1825  const dof_id_type their_dof_g = neigh_dof_indices[i];
1826  libmesh_assert_not_equal_to (their_dof_g, DofObject::invalid_id);
1827  libmesh_assert_not_equal_to (their_dof_g, my_dof_g);
1828 
1829  const Real their_dof_value = Ue[is](js);
1830 
1831  if (std::abs(their_dof_value) < 10*TOLERANCE)
1832  continue;
1833 
1834  constraint_row->insert(std::make_pair(their_dof_g,
1835  their_dof_value));
1836  }
1837  }
1838  }
1839  // p refinement constraints:
1840  // constrain dofs shared between
1841  // active elements and neighbors with
1842  // lower polynomial degrees
1843  const unsigned int min_p_level =
1844  neigh->min_p_level_by_neighbor(elem, elem->p_level());
1845  if (min_p_level < elem->p_level())
1846  {
1847  // Adaptive p refinement of non-hierarchic bases will
1848  // require more coding
1849  libmesh_assert(my_fe->is_hierarchic());
1850  dof_map.constrain_p_dofs(variable_number, elem,
1851  s, min_p_level);
1852  }
1853  }
1854 }
1855 
1856 #endif // #ifdef LIBMESH_ENABLE_AMR
1857 
1858 
1859 
1860 #ifdef LIBMESH_ENABLE_PERIODIC
1861 template <typename OutputType>
1862 void
1865  DofMap &dof_map,
1866  const PeriodicBoundaries &boundaries,
1867  const MeshBase &mesh,
1868  const PointLocatorBase *point_locator,
1869  const unsigned int variable_number,
1870  const Elem* elem)
1871 {
1872  // Only bother if we truly have periodic boundaries
1873  if (boundaries.empty())
1874  return;
1875 
1876  libmesh_assert(elem);
1877 
1878  // Only constrain active elements with this method
1879  if (!elem->active())
1880  return;
1881 
1882  const unsigned int Dim = elem->dim();
1883 
1884  // We need sys_number and variable_number for DofObject methods
1885  // later
1886  const unsigned int sys_number = dof_map.sys_number();
1887 
1888  const FEType& base_fe_type = dof_map.variable_type(variable_number);
1889 
1890  // Construct FE objects for this element and its pseudo-neighbors.
1892  (FEGenericBase<OutputShape>::build(Dim, base_fe_type));
1893  const FEContinuity cont = my_fe->get_continuity();
1894 
1895  // We don't need to constrain discontinuous elements
1896  if (cont == DISCONTINUOUS)
1897  return;
1898  libmesh_assert (cont == C_ZERO || cont == C_ONE);
1899 
1900  // We'll use element size to generate relative tolerances later
1901  const Real primary_hmin = elem->hmin();
1902 
1904  (FEGenericBase<OutputShape>::build(Dim, base_fe_type));
1905 
1906  QGauss my_qface(Dim-1, base_fe_type.default_quadrature_order());
1907  my_fe->attach_quadrature_rule (&my_qface);
1908  std::vector<Point> neigh_qface;
1909 
1910  const std::vector<Real>& JxW = my_fe->get_JxW();
1911  const std::vector<Point>& q_point = my_fe->get_xyz();
1912  const std::vector<std::vector<OutputShape> >& phi = my_fe->get_phi();
1913  const std::vector<std::vector<OutputShape> >& neigh_phi =
1914  neigh_fe->get_phi();
1915  const std::vector<Point> *face_normals = NULL;
1916  const std::vector<std::vector<OutputGradient> > *dphi = NULL;
1917  const std::vector<std::vector<OutputGradient> > *neigh_dphi = NULL;
1918  std::vector<dof_id_type> my_dof_indices, neigh_dof_indices;
1919  std::vector<unsigned int> my_side_dofs, neigh_side_dofs;
1920 
1921  if (cont != C_ZERO)
1922  {
1923  const std::vector<Point>& ref_face_normals =
1924  my_fe->get_normals();
1925  face_normals = &ref_face_normals;
1926  const std::vector<std::vector<OutputGradient> >& ref_dphi =
1927  my_fe->get_dphi();
1928  dphi = &ref_dphi;
1929  const std::vector<std::vector<OutputGradient> >& ref_neigh_dphi =
1930  neigh_fe->get_dphi();
1931  neigh_dphi = &ref_neigh_dphi;
1932  }
1933 
1934  DenseMatrix<Real> Ke;
1935  DenseVector<Real> Fe;
1936  std::vector<DenseVector<Real> > Ue;
1937 
1938  // Look at the element faces. Check to see if we need to
1939  // build constraints.
1940  for (unsigned int s=0; s<elem->n_sides(); s++)
1941  {
1942  if (elem->neighbor(s))
1943  continue;
1944 
1945  const std::vector<boundary_id_type>& bc_ids = mesh.boundary_info->boundary_ids (elem, s);
1946  for (std::vector<boundary_id_type>::const_iterator id_it=bc_ids.begin(); id_it!=bc_ids.end(); ++id_it)
1947  {
1948  const boundary_id_type boundary_id = *id_it;
1949  const PeriodicBoundaryBase *periodic = boundaries.boundary(boundary_id);
1950  if (periodic && periodic->is_my_variable(variable_number))
1951  {
1952  libmesh_assert(point_locator);
1953 
1954  // Get pointers to the element's neighbor.
1955  const Elem* neigh = boundaries.neighbor(boundary_id, *point_locator, elem, s);
1956 
1957  if (neigh == NULL)
1958  {
1959  libMesh::err << "PeriodicBoundaries point locator object returned NULL!" << std::endl;
1960  libmesh_error();
1961  }
1962 
1963  // periodic (and possibly h refinement) constraints:
1964  // constrain dofs shared between
1965  // this element and ones as coarse
1966  // as or coarser than this element.
1967  if (neigh->level() <= elem->level())
1968  {
1969  unsigned int s_neigh =
1970  mesh.boundary_info->side_with_boundary_id (neigh, periodic->pairedboundary);
1971  libmesh_assert_not_equal_to (s_neigh, libMesh::invalid_uint);
1972 
1973 #ifdef LIBMESH_ENABLE_AMR
1974  // Find the minimum p level; we build the h constraint
1975  // matrix with this and then constrain away all higher p
1976  // DoFs.
1977  libmesh_assert(neigh->active());
1978  const unsigned int min_p_level =
1979  std::min(elem->p_level(), neigh->p_level());
1980 
1981  // we may need to make the FE objects reinit with the
1982  // minimum shared p_level
1983  // FIXME - I hate using const_cast<> and avoiding
1984  // accessor functions; there's got to be a
1985  // better way to do this!
1986  const unsigned int old_elem_level = elem->p_level();
1987  if (old_elem_level != min_p_level)
1988  (const_cast<Elem *>(elem))->hack_p_level(min_p_level);
1989  const unsigned int old_neigh_level = neigh->p_level();
1990  if (old_neigh_level != min_p_level)
1991  (const_cast<Elem *>(neigh))->hack_p_level(min_p_level);
1992 #endif // #ifdef LIBMESH_ENABLE_AMR
1993 
1994  // We can do a projection with a single integration,
1995  // due to the assumption of nested finite element
1996  // subspaces.
1997  // FIXME: it might be more efficient to do nodes,
1998  // then edges, then side, to reduce the size of the
1999  // Cholesky factorization(s)
2000  my_fe->reinit(elem, s);
2001 
2002  dof_map.dof_indices (elem, my_dof_indices,
2003  variable_number);
2004  dof_map.dof_indices (neigh, neigh_dof_indices,
2005  variable_number);
2006 
2007  const unsigned int n_qp = my_qface.n_points();
2008 
2009  // Translate the quadrature points over to the
2010  // neighbor's boundary
2011  std::vector<Point> neigh_point(q_point.size());
2012  for (unsigned int i=0; i != neigh_point.size(); ++i)
2013  neigh_point[i] = periodic->get_corresponding_pos(q_point[i]);
2014 
2015  FEInterface::inverse_map (Dim, base_fe_type, neigh,
2016  neigh_point, neigh_qface);
2017 
2018  neigh_fe->reinit(neigh, &neigh_qface);
2019 
2020  // We're only concerned with DOFs whose values (and/or first
2021  // derivatives for C1 elements) are supported on side nodes
2022  FEInterface::dofs_on_side(elem, Dim, base_fe_type, s, my_side_dofs);
2023  FEInterface::dofs_on_side(neigh, Dim, base_fe_type, s_neigh, neigh_side_dofs);
2024 
2025  // We're done with functions that examine Elem::p_level(),
2026  // so let's unhack those levels
2027 #ifdef LIBMESH_ENABLE_AMR
2028  if (elem->p_level() != old_elem_level)
2029  (const_cast<Elem *>(elem))->hack_p_level(old_elem_level);
2030  if (neigh->p_level() != old_neigh_level)
2031  (const_cast<Elem *>(neigh))->hack_p_level(old_neigh_level);
2032 #endif // #ifdef LIBMESH_ENABLE_AMR
2033 
2034  const unsigned int n_side_dofs =
2035  libmesh_cast_int<unsigned int>
2036  (my_side_dofs.size());
2037  libmesh_assert_equal_to (n_side_dofs, neigh_side_dofs.size());
2038 
2039  Ke.resize (n_side_dofs, n_side_dofs);
2040  Ue.resize(n_side_dofs);
2041 
2042  // Form the projection matrix, (inner product of fine basis
2043  // functions against fine test functions)
2044  for (unsigned int is = 0; is != n_side_dofs; ++is)
2045  {
2046  const unsigned int i = my_side_dofs[is];
2047  for (unsigned int js = 0; js != n_side_dofs; ++js)
2048  {
2049  const unsigned int j = my_side_dofs[js];
2050  for (unsigned int qp = 0; qp != n_qp; ++qp)
2051  {
2052  Ke(is,js) += JxW[qp] *
2053  TensorTools::inner_product(phi[i][qp],
2054  phi[j][qp]);
2055  if (cont != C_ZERO)
2056  Ke(is,js) += JxW[qp] *
2057  TensorTools::inner_product((*dphi)[i][qp] *
2058  (*face_normals)[qp],
2059  (*dphi)[j][qp] *
2060  (*face_normals)[qp]);
2061  }
2062  }
2063  }
2064 
2065  // Form the right hand sides, (inner product of coarse basis
2066  // functions against fine test functions)
2067  for (unsigned int is = 0; is != n_side_dofs; ++is)
2068  {
2069  const unsigned int i = neigh_side_dofs[is];
2070  Fe.resize (n_side_dofs);
2071  for (unsigned int js = 0; js != n_side_dofs; ++js)
2072  {
2073  const unsigned int j = my_side_dofs[js];
2074  for (unsigned int qp = 0; qp != n_qp; ++qp)
2075  {
2076  Fe(js) += JxW[qp] *
2077  TensorTools::inner_product(neigh_phi[i][qp],
2078  phi[j][qp]);
2079  if (cont != C_ZERO)
2080  Fe(js) += JxW[qp] *
2081  TensorTools::inner_product((*neigh_dphi)[i][qp] *
2082  (*face_normals)[qp],
2083  (*dphi)[j][qp] *
2084  (*face_normals)[qp]);
2085  }
2086  }
2087  Ke.cholesky_solve(Fe, Ue[is]);
2088  }
2089 
2090  // Make sure we're not adding recursive constraints
2091  // due to the redundancy in the way we add periodic
2092  // boundary constraints
2093  //
2094  // In order for this to work while threaded or on
2095  // distributed meshes, we need a rigorous way to
2096  // avoid recursive constraints. Here it is:
2097  //
2098  // For vertex DoFs, if there is a "prior" element
2099  // (i.e. a coarser element or an equally refined
2100  // element with a lower id) on this boundary which
2101  // contains the vertex point, then we will avoid
2102  // generating constraints; the prior element (or
2103  // something prior to it) may do so. If we are the
2104  // most prior (or "primary") element on this
2105  // boundary sharing this point, then we look at the
2106  // boundary periodic to us, we find the primary
2107  // element there, and if that primary is coarser or
2108  // equal-but-lower-id, then our vertex dofs are
2109  // constrained in terms of that element.
2110  //
2111  // For edge DoFs, if there is a coarser element
2112  // on this boundary sharing this edge, then we will
2113  // avoid generating constraints (we will be
2114  // constrained indirectly via AMR constraints
2115  // connecting us to the coarser element's DoFs). If
2116  // we are the coarsest element sharing this edge,
2117  // then we generate constraints if and only if we
2118  // are finer than the coarsest element on the
2119  // boundary periodic to us sharing the corresponding
2120  // periodic edge, or if we are at equal level but
2121  // our edge nodes have higher ids than the periodic
2122  // edge nodes (sorted from highest to lowest, then
2123  // compared lexicographically)
2124  //
2125  // For face DoFs, we generate constraints if we are
2126  // finer than our periodic neighbor, or if we are at
2127  // equal level but our element id is higher than its
2128  // element id.
2129  //
2130  // If the primary neighbor is also the current elem
2131  // (a 1-element-thick mesh) then we choose which
2132  // vertex dofs to constrain via lexicographic
2133  // ordering on point locations
2134 
2135  // FIXME: This code doesn't yet properly handle
2136  // cases where multiple different periodic BCs
2137  // intersect.
2138  std::set<dof_id_type> my_constrained_dofs;
2139 
2140  for (unsigned int n = 0; n != elem->n_nodes(); ++n)
2141  {
2142  if (!elem->is_node_on_side(n,s))
2143  continue;
2144 
2145  const Node* my_node = elem->get_node(n);
2146 
2147  if (elem->is_vertex(n))
2148  {
2149  // Find all boundary ids that include this
2150  // point and have periodic boundary
2151  // conditions for this variable
2152  std::set<boundary_id_type> point_bcids;
2153 
2154  for (unsigned int new_s = 0; new_s !=
2155  elem->n_sides(); ++new_s)
2156  {
2157  if (!elem->is_node_on_side(n,new_s))
2158  continue;
2159 
2160  const std::vector<boundary_id_type>
2161  new_bc_ids = mesh.boundary_info->boundary_ids (elem, s);
2162  for (std::vector<boundary_id_type>::const_iterator
2163  new_id_it=new_bc_ids.begin(); new_id_it!=new_bc_ids.end(); ++new_id_it)
2164  {
2165  const boundary_id_type new_boundary_id = *new_id_it;
2166  const PeriodicBoundaryBase *new_periodic = boundaries.boundary(new_boundary_id);
2167  if (new_periodic && new_periodic->is_my_variable(variable_number))
2168  {
2169  point_bcids.insert(new_boundary_id);
2170  }
2171  }
2172  }
2173 
2174  // See if this vertex has point neighbors to
2175  // defer to
2176  if (primary_boundary_point_neighbor
2177  (elem, *my_node, *mesh.boundary_info, point_bcids) != elem)
2178  continue;
2179 
2180  // Find the complementary boundary id set
2181  std::set<boundary_id_type> point_pairedids;
2182  for (std::set<boundary_id_type>::const_iterator i =
2183  point_bcids.begin(); i != point_bcids.end(); ++i)
2184  {
2185  const boundary_id_type new_boundary_id = *i;
2186  const PeriodicBoundaryBase *new_periodic = boundaries.boundary(new_boundary_id);
2187  point_pairedids.insert(new_periodic->pairedboundary);
2188  }
2189 
2190  // What do we want to constrain against?
2191  const Elem* primary_elem = NULL;
2192  const Elem* main_neigh = NULL;
2193  Point main_pt = *my_node,
2194  primary_pt = *my_node;
2195 
2196  for (std::set<boundary_id_type>::const_iterator i =
2197  point_bcids.begin(); i != point_bcids.end(); ++i)
2198  {
2199  // Find the corresponding periodic point and
2200  // its primary neighbor
2201  const boundary_id_type new_boundary_id = *i;
2202  const PeriodicBoundaryBase *new_periodic = boundaries.boundary(new_boundary_id);
2203 
2204  const Point neigh_pt =
2205  new_periodic->get_corresponding_pos(*my_node);
2206 
2207  // If the point is getting constrained
2208  // to itself by this PBC then we don't
2209  // generate any constraints
2210  if (neigh_pt.absolute_fuzzy_equals
2211  (*my_node, primary_hmin*TOLERANCE))
2212  continue;
2213 
2214  // Otherwise we'll have a constraint in
2215  // one direction or another
2216  if (!primary_elem)
2217  primary_elem = elem;
2218 
2219  const Elem *primary_neigh = primary_boundary_point_neighbor
2220  (neigh, neigh_pt, *mesh.boundary_info,
2221  point_pairedids);
2222 
2223  libmesh_assert(primary_neigh);
2224 
2225  if (new_boundary_id == boundary_id)
2226  {
2227  main_neigh = primary_neigh;
2228  main_pt = neigh_pt;
2229  }
2230 
2231  // Finer elements will get constrained in
2232  // terms of coarser neighbors, not the
2233  // other way around
2234  if ((primary_neigh->level() > primary_elem->level()) ||
2235 
2236  // For equal-level elements, the one with
2237  // higher id gets constrained in terms of
2238  // the one with lower id
2239  (primary_neigh->level() == primary_elem->level() &&
2240  primary_neigh->id() > primary_elem->id()) ||
2241 
2242  // On a one-element-thick mesh, we compare
2243  // points to see what side gets constrained
2244  (primary_neigh == primary_elem &&
2245  (neigh_pt > primary_pt)))
2246  continue;
2247 
2248  primary_elem = primary_neigh;
2249  primary_pt = neigh_pt;
2250  }
2251 
2252  if (!primary_elem ||
2253  primary_elem != main_neigh ||
2254  primary_pt != main_pt)
2255  continue;
2256  }
2257  else if (elem->is_edge(n))
2258  {
2259  // Find which edge we're on
2260  unsigned int e=0;
2261  for (; e != elem->n_edges(); ++e)
2262  {
2263  if (elem->is_node_on_edge(n,e))
2264  break;
2265  }
2266  libmesh_assert_less (e, elem->n_edges());
2267 
2268  // Find the edge end nodes
2269  Node *e1 = NULL,
2270  *e2 = NULL;
2271  for (unsigned int nn = 0; nn != elem->n_nodes(); ++nn)
2272  {
2273  if (nn == n)
2274  continue;
2275 
2276  if (elem->is_node_on_edge(nn, e))
2277  {
2278  if (e1 == NULL)
2279  {
2280  e1 = elem->get_node(nn);
2281  }
2282  else
2283  {
2284  e2 = elem->get_node(nn);
2285  break;
2286  }
2287  }
2288  }
2289  libmesh_assert (e1 && e2);
2290 
2291  // Find all boundary ids that include this
2292  // edge and have periodic boundary
2293  // conditions for this variable
2294  std::set<boundary_id_type> edge_bcids;
2295 
2296  for (unsigned int new_s = 0; new_s !=
2297  elem->n_sides(); ++new_s)
2298  {
2299  if (!elem->is_node_on_side(n,new_s))
2300  continue;
2301 
2302  const std::vector<boundary_id_type>&
2303  new_bc_ids = mesh.boundary_info->boundary_ids (elem, s);
2304  for (std::vector<boundary_id_type>::const_iterator
2305  new_id_it=new_bc_ids.begin(); new_id_it!=new_bc_ids.end(); ++new_id_it)
2306  {
2307  const boundary_id_type new_boundary_id = *new_id_it;
2308  const PeriodicBoundaryBase *new_periodic = boundaries.boundary(new_boundary_id);
2309  if (new_periodic && new_periodic->is_my_variable(variable_number))
2310  {
2311  edge_bcids.insert(new_boundary_id);
2312  }
2313  }
2314  }
2315 
2316 
2317  // See if this edge has neighbors to defer to
2318  if (primary_boundary_edge_neighbor
2319  (elem, *e1, *e2, *mesh.boundary_info, edge_bcids) != elem)
2320  continue;
2321 
2322  // Find the complementary boundary id set
2323  std::set<boundary_id_type> edge_pairedids;
2324  for (std::set<boundary_id_type>::const_iterator i =
2325  edge_bcids.begin(); i != edge_bcids.end(); ++i)
2326  {
2327  const boundary_id_type new_boundary_id = *i;
2328  const PeriodicBoundaryBase *new_periodic = boundaries.boundary(new_boundary_id);
2329  edge_pairedids.insert(new_periodic->pairedboundary);
2330  }
2331 
2332  // What do we want to constrain against?
2333  const Elem* primary_elem = NULL;
2334  const Elem* main_neigh = NULL;
2335  Point main_pt1 = *e1,
2336  main_pt2 = *e2,
2337  primary_pt1 = *e1,
2338  primary_pt2 = *e2;
2339 
2340  for (std::set<boundary_id_type>::const_iterator i =
2341  edge_bcids.begin(); i != edge_bcids.end(); ++i)
2342  {
2343  // Find the corresponding periodic edge and
2344  // its primary neighbor
2345  const boundary_id_type new_boundary_id = *i;
2346  const PeriodicBoundaryBase *new_periodic = boundaries.boundary(new_boundary_id);
2347 
2348  Point neigh_pt1 = new_periodic->get_corresponding_pos(*e1),
2349  neigh_pt2 = new_periodic->get_corresponding_pos(*e2);
2350 
2351  // If the edge is getting constrained
2352  // to itself by this PBC then we don't
2353  // generate any constraints
2354  if (neigh_pt1.absolute_fuzzy_equals
2355  (*e1, primary_hmin*TOLERANCE) &&
2356  neigh_pt2.absolute_fuzzy_equals
2357  (*e2, primary_hmin*TOLERANCE))
2358  continue;
2359 
2360  // Otherwise we'll have a constraint in
2361  // one direction or another
2362  if (!primary_elem)
2363  primary_elem = elem;
2364 
2365  const Elem *primary_neigh = primary_boundary_edge_neighbor
2366  (neigh, neigh_pt1, neigh_pt2, *mesh.boundary_info,
2367  edge_pairedids);
2368 
2369  libmesh_assert(primary_neigh);
2370 
2371  if (new_boundary_id == boundary_id)
2372  {
2373  main_neigh = primary_neigh;
2374  main_pt1 = neigh_pt1;
2375  main_pt2 = neigh_pt2;
2376  }
2377 
2378  // If we have a one-element thick mesh,
2379  // we'll need to sort our points to get a
2380  // consistent ordering rule
2381  //
2382  // Use >= in this test to make sure that,
2383  // for angular constraints, no node gets
2384  // constrained to itself.
2385  if (primary_neigh == primary_elem)
2386  {
2387  if (primary_pt1 > primary_pt2)
2388  std::swap(primary_pt1, primary_pt2);
2389  if (neigh_pt1 > neigh_pt2)
2390  std::swap(neigh_pt1, neigh_pt2);
2391 
2392  if (neigh_pt2 >= primary_pt2)
2393  continue;
2394  }
2395 
2396  // Otherwise:
2397  // Finer elements will get constrained in
2398  // terms of coarser ones, not the other way
2399  // around
2400  if ((primary_neigh->level() > primary_elem->level()) ||
2401 
2402  // For equal-level elements, the one with
2403  // higher id gets constrained in terms of
2404  // the one with lower id
2405  (primary_neigh->level() == primary_elem->level() &&
2406  primary_neigh->id() > primary_elem->id()))
2407  continue;
2408 
2409  primary_elem = primary_neigh;
2410  primary_pt1 = neigh_pt1;
2411  primary_pt2 = neigh_pt2;
2412  }
2413 
2414  if (!primary_elem ||
2415  primary_elem != main_neigh ||
2416  primary_pt1 != main_pt1 ||
2417  primary_pt2 != main_pt2)
2418  continue;
2419  }
2420  else if (elem->is_face(n))
2421  {
2422  // If we have a one-element thick mesh,
2423  // use the ordering of the face node and its
2424  // periodic counterpart to determine what
2425  // gets constrained
2426  if (neigh == elem)
2427  {
2428  const Point neigh_pt =
2429  periodic->get_corresponding_pos(*my_node);
2430  if (neigh_pt > *my_node)
2431  continue;
2432  }
2433 
2434  // Otherwise:
2435  // Finer elements will get constrained in
2436  // terms of coarser ones, not the other way
2437  // around
2438  if ((neigh->level() > elem->level()) ||
2439 
2440  // For equal-level elements, the one with
2441  // higher id gets constrained in terms of
2442  // the one with lower id
2443  (neigh->level() == elem->level() &&
2444  neigh->id() > elem->id()))
2445  continue;
2446  }
2447 
2448  // If we made it here without hitting a continue
2449  // statement, then we're at a node whose dofs
2450  // should be constrained by this element's
2451  // calculations.
2452  const unsigned int n_comp =
2453  my_node->n_comp(sys_number, variable_number);
2454 
2455  for (unsigned int i=0; i != n_comp; ++i)
2456  my_constrained_dofs.insert
2457  (my_node->dof_number
2458  (sys_number, variable_number, i));
2459  }
2460 
2461  // FIXME: old code for disambiguating periodic BCs:
2462  // this is not threadsafe nor safe to run on a
2463  // non-serialized mesh.
2464  /*
2465  std::vector<bool> recursive_constraint(n_side_dofs, false);
2466 
2467  for (unsigned int is = 0; is != n_side_dofs; ++is)
2468  {
2469  const unsigned int i = neigh_side_dofs[is];
2470  const dof_id_type their_dof_g = neigh_dof_indices[i];
2471  libmesh_assert_not_equal_to (their_dof_g, DofObject::invalid_id);
2472 
2473  {
2474  Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
2475 
2476  if (!dof_map.is_constrained_dof(their_dof_g))
2477  continue;
2478  }
2479 
2480  DofConstraintRow& their_constraint_row =
2481  constraints[their_dof_g].first;
2482 
2483  for (unsigned int js = 0; js != n_side_dofs; ++js)
2484  {
2485  const unsigned int j = my_side_dofs[js];
2486  const dof_id_type my_dof_g = my_dof_indices[j];
2487  libmesh_assert_not_equal_to (my_dof_g, DofObject::invalid_id);
2488 
2489  if (their_constraint_row.count(my_dof_g))
2490  recursive_constraint[js] = true;
2491  }
2492  }
2493  */
2494 
2495  for (unsigned int js = 0; js != n_side_dofs; ++js)
2496  {
2497  // FIXME: old code path
2498  // if (recursive_constraint[js])
2499  // continue;
2500 
2501  const unsigned int j = my_side_dofs[js];
2502  const dof_id_type my_dof_g = my_dof_indices[j];
2503  libmesh_assert_not_equal_to (my_dof_g, DofObject::invalid_id);
2504 
2505  // FIXME: new code path
2506  if (!my_constrained_dofs.count(my_dof_g))
2507  continue;
2508 
2509  DofConstraintRow* constraint_row;
2510 
2511  // we may be running constraint methods concurretly
2512  // on multiple threads, so we need a lock to
2513  // ensure that this constraint is "ours"
2514  {
2515  Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
2516 
2517  if (dof_map.is_constrained_dof(my_dof_g))
2518  continue;
2519 
2520  constraint_row = &(constraints[my_dof_g]);
2521  libmesh_assert(constraint_row->empty());
2522  }
2523 
2524  for (unsigned int is = 0; is != n_side_dofs; ++is)
2525  {
2526  const unsigned int i = neigh_side_dofs[is];
2527  const dof_id_type their_dof_g = neigh_dof_indices[i];
2528  libmesh_assert_not_equal_to (their_dof_g, DofObject::invalid_id);
2529 
2530  // Periodic constraints should never be
2531  // self-constraints
2532  // libmesh_assert_not_equal_to (their_dof_g, my_dof_g);
2533 
2534  const Real their_dof_value = Ue[is](js);
2535 
2536  if (their_dof_g == my_dof_g)
2537  {
2538  libmesh_assert_less (std::abs(their_dof_value-1.), 1.e-5);
2539  for (unsigned int k = 0; k != n_side_dofs; ++k)
2540  libmesh_assert(k == is || std::abs(Ue[k](js)) < 1.e-5);
2541  continue;
2542  }
2543 
2544  if (std::abs(their_dof_value) < 10*TOLERANCE)
2545  continue;
2546 
2547  constraint_row->insert(std::make_pair(their_dof_g,
2548  their_dof_value));
2549  }
2550  }
2551  }
2552  // p refinement constraints:
2553  // constrain dofs shared between
2554  // active elements and neighbors with
2555  // lower polynomial degrees
2556 #ifdef LIBMESH_ENABLE_AMR
2557  const unsigned int min_p_level =
2558  neigh->min_p_level_by_neighbor(elem, elem->p_level());
2559  if (min_p_level < elem->p_level())
2560  {
2561  // Adaptive p refinement of non-hierarchic bases will
2562  // require more coding
2563  libmesh_assert(my_fe->is_hierarchic());
2564  dof_map.constrain_p_dofs(variable_number, elem,
2565  s, min_p_level);
2566  }
2567 #endif // #ifdef LIBMESH_ENABLE_AMR
2568  }
2569  }
2570  }
2571 }
2572 
2573 #endif // LIBMESH_ENABLE_PERIODIC
2574 
2575 // ------------------------------------------------------------
2576 // Explicit instantiations
2577 template class FEGenericBase<Real>;
2578 template class FEGenericBase<RealGradient>;
2579 
2580 } // namespace libMesh

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

Hosted By:
SourceForge.net Logo