00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012 #ifndef __deal2__fe_h
00013 #define __deal2__fe_h
00014
00015 #include <deal.II/base/config.h>
00016 #include <deal.II/base/geometry_info.h>
00017 #include <deal.II/fe/fe_base.h>
00018
00019
00020 DEAL_II_NAMESPACE_OPEN
00021
00022 template <int dim, int spacedim> class FEValuesData;
00023 template <int dim, int spacedim> class FEValuesBase;
00024 template <int dim, int spacedim> class FEValues;
00025 template <int dim, int spacedim> class FEFaceValues;
00026 template <int dim, int spacedim> class FESubfaceValues;
00027 template <int dim, int spacedim> class FESystem;
00028 template <class POLY, int dim, int spacedim> class FE_PolyTensor;
00029 namespace hp
00030 {
00031 template <int dim, int spacedim> class FECollection;
00032 }
00033
00034
00335 template <int dim, int spacedim=dim>
00336 class FiniteElement : public Subscriptor,
00337 public FiniteElementData<dim>
00338 {
00339 public:
00352 class InternalDataBase : public Mapping<dim,spacedim>::InternalDataBase
00353 {
00354 public:
00360 virtual ~InternalDataBase ();
00361
00369 void initialize_2nd (const FiniteElement<dim,spacedim> *element,
00370 const Mapping<dim,spacedim> &mapping,
00371 const Quadrature<dim> &quadrature);
00372
00386 std::vector<FEValues<dim,spacedim>*> differences;
00387 };
00388
00389 public:
00393 FiniteElement (const FiniteElementData<dim> &fe_data,
00394 const std::vector<bool> &restriction_is_additive_flags,
00395 const std::vector<std::vector<bool> > &nonzero_components);
00396
00402 virtual ~FiniteElement ();
00403
00422 virtual std::string get_name () const = 0;
00423
00467 const FiniteElement<dim,spacedim> & operator[] (const unsigned int fe_index) const;
00468
00504 virtual double shape_value (const unsigned int i,
00505 const Point<dim> &p) const;
00506
00519 virtual double shape_value_component (const unsigned int i,
00520 const Point<dim> &p,
00521 const unsigned int component) const;
00522
00557 virtual Tensor<1,dim> shape_grad (const unsigned int i,
00558 const Point<dim> &p) const;
00559
00572 virtual Tensor<1,dim> shape_grad_component (const unsigned int i,
00573 const Point<dim> &p,
00574 const unsigned int component) const;
00575
00609 virtual Tensor<2,dim> shape_grad_grad (const unsigned int i,
00610 const Point<dim> &p) const;
00611
00624 virtual Tensor<2,dim> shape_grad_grad_component (const unsigned int i,
00625 const Point<dim> &p,
00626 const unsigned int component) const;
00644 virtual bool has_support_on_face (const unsigned int shape_index,
00645 const unsigned int face_index) const;
00646
00648
00692 const FullMatrix<double> &
00693 get_restriction_matrix (const unsigned int child,
00694 const RefinementCase<dim> &refinement_case=RefinementCase<dim>::isotropic_refinement) const;
00695
00748 const FullMatrix<double> &
00749 get_prolongation_matrix (const unsigned int child,
00750 const RefinementCase<dim> &refinement_case=RefinementCase<dim>::isotropic_refinement) const;
00751
00790 bool prolongation_is_implemented () const;
00791
00818 bool isotropic_prolongation_is_implemented () const;
00819
00858 bool restriction_is_implemented () const;
00859
00886 bool isotropic_restriction_is_implemented () const;
00887
00888
00899 bool restriction_is_additive (const unsigned int index) const;
00900
00929 const FullMatrix<double> & constraints (const internal::SubfaceCase<dim> &subface_case=internal::SubfaceCase<dim>::case_isotropic) const;
00930
00961 bool constraints_are_implemented (const internal::SubfaceCase<dim> &subface_case=internal::SubfaceCase<dim>::case_isotropic) const;
00962
00963
01001 virtual bool hp_constraints_are_implemented () const;
01002
01003
01024 virtual void
01025 get_interpolation_matrix (const FiniteElement<dim,spacedim> &source,
01026 FullMatrix<double> &matrix) const;
01028
01056 virtual void
01057 get_face_interpolation_matrix (const FiniteElement<dim,spacedim> &source,
01058 FullMatrix<double> &matrix) const;
01059
01060
01082 virtual void
01083 get_subface_interpolation_matrix (const FiniteElement<dim,spacedim> &source,
01084 const unsigned int subface,
01085 FullMatrix<double> &matrix) const;
01087
01088
01123 virtual
01124 std::vector<std::pair<unsigned int, unsigned int> >
01125 hp_vertex_dof_identities (const FiniteElement<dim,spacedim> &fe_other) const;
01126
01134 virtual
01135 std::vector<std::pair<unsigned int, unsigned int> >
01136 hp_line_dof_identities (const FiniteElement<dim,spacedim> &fe_other) const;
01137
01145 virtual
01146 std::vector<std::pair<unsigned int, unsigned int> >
01147 hp_quad_dof_identities (const FiniteElement<dim,spacedim> &fe_other) const;
01148
01161 virtual
01162 FiniteElementDomination::Domination
01163 compare_for_face_domination (const FiniteElement<dim,spacedim> &fe_other) const;
01164
01166
01184 bool operator == (const FiniteElement<dim,spacedim> &) const;
01185
01247 std::pair<unsigned int, unsigned int>
01248 system_to_component_index (const unsigned int index) const;
01249
01265 unsigned int component_to_system_index(const unsigned int component,
01266 const unsigned int index) const;
01267
01286 std::pair<unsigned int, unsigned int>
01287 face_system_to_component_index (const unsigned int index) const;
01288
01302 unsigned int adjust_quad_dof_index_for_face_orientation (const unsigned int index,
01303 const bool face_orientation,
01304 const bool face_flip,
01305 const bool face_rotation) const;
01306
01319 unsigned int adjust_line_dof_index_for_line_orientation (const unsigned int index,
01320 const bool line_orientation) const;
01321
01352 const std::vector<bool> &
01353 get_nonzero_components (const unsigned int i) const;
01354
01375 unsigned int
01376 n_nonzero_components (const unsigned int i) const;
01377
01397 bool
01398 is_primitive (const unsigned int i) const;
01399
01405 using FiniteElementData<dim>::is_primitive;
01406
01429 unsigned int n_base_elements () const;
01430
01438 virtual
01439 const FiniteElement<dim,spacedim> &
01440 base_element (const unsigned int index) const;
01441
01452 unsigned int
01453 element_multiplicity (const unsigned int index) const;
01454
01492 std::pair<std::pair<unsigned int, unsigned int>, unsigned int>
01493 system_to_base_index (const unsigned int index) const;
01494
01512 std::pair<std::pair<unsigned int, unsigned int>, unsigned int>
01513 face_system_to_base_index (const unsigned int index) const;
01514
01520 unsigned int first_block_of_base (const unsigned int b) const;
01521
01546 std::pair<unsigned int, unsigned int>
01547 component_to_base_index (const unsigned int component) const;
01548
01549
01555 std::pair<unsigned int,unsigned int>
01556 block_to_base_index (const unsigned int block) const;
01557
01563 std::pair<unsigned int,unsigned int>
01564 system_to_block_index (const unsigned int component) const;
01565
01570 unsigned int
01571 component_to_block_index (const unsigned int component) const;
01572
01574
01609 const std::vector<Point<dim> > &
01610 get_unit_support_points () const;
01611
01639 bool has_support_points () const;
01640
01667 virtual
01668 Point<dim>
01669 unit_support_point (const unsigned int index) const;
01670
01727 const std::vector<Point<dim-1> > &
01728 get_unit_face_support_points () const;
01729
01743 bool has_face_support_points () const;
01744
01751 virtual
01752 Point<dim-1>
01753 unit_face_support_point (const unsigned int index) const;
01754
01759 const std::vector<Point<dim> > &
01760 get_generalized_support_points () const;
01761
01769 bool has_generalized_support_points () const;
01770
01774 const std::vector<Point<dim-1> > &
01775 get_generalized_face_support_points () const;
01776
01791 bool has_generalized_face_support_points () const;
01792
01808 virtual void interpolate(std::vector<double>& local_dofs,
01809 const std::vector<double>& values) const;
01810
01825 virtual void interpolate(std::vector<double>& local_dofs,
01826 const std::vector<Vector<double> >& values,
01827 unsigned int offset = 0) const;
01828
01834 virtual void interpolate(
01835 std::vector<double>& local_dofs,
01836 const VectorSlice<const std::vector<std::vector<double> > >& values) const;
01837
01839
01851 virtual std::size_t memory_consumption () const;
01857 DeclException1 (ExcShapeFunctionNotPrimitive,
01858 int,
01859 << "The shape function with index " << arg1
01860 << " is not primitive, i.e. it is vector-valued and "
01861 << "has more than one non-zero vector component. This "
01862 << "function cannot be called for these shape functions. "
01863 << "Maybe you want to use the same function with the "
01864 << "_component suffix?");
01870 DeclException0 (ExcFENotPrimitive);
01876 DeclException0 (ExcUnitShapeValuesDoNotExist);
01877
01885 DeclException0 (ExcFEHasNoSupportPoints);
01886
01895 DeclException0 (ExcEmbeddingVoid);
01896
01906 DeclException0 (ExcProjectionVoid);
01907
01917 DeclException0 (ExcConstraintsVoid);
01918
01923 DeclException2 (ExcWrongInterfaceMatrixSize,
01924 int, int,
01925 << "The interface matrix has a size of " << arg1
01926 << "x" << arg2
01927 << ", which is not reasonable in the present dimension.");
01932 DeclException2 (ExcComponentIndexInvalid,
01933 int, int,
01934 << "The component-index pair (" << arg1 << ", " << arg2
01935 << ") is invalid, i.e. non-existent");
01940 DeclException0 (ExcInterpolationNotImplemented);
01941
01947 DeclException0 (ExcBoundaryFaceUsed);
01953 DeclException0 (ExcJacobiDeterminantHasWrongSign);
01954
01955 protected:
01956
01987 void reinit_restriction_and_prolongation_matrices(const bool isotropic_restriction_only=false,
01988 const bool isotropic_prolongation_only=false);
01989
02016 std::vector<std::vector<FullMatrix<double> > > restriction;
02017
02044 std::vector<std::vector<FullMatrix<double> > > prolongation;
02045
02061 FullMatrix<double> interface_constraints;
02062
02082 std::vector<Point<dim> > unit_support_points;
02083
02092 std::vector<Point<dim-1> > unit_face_support_points;
02093
02099 std::vector<Point<dim> > generalized_support_points;
02100
02106 std::vector<Point<dim-1> > generalized_face_support_points;
02107
02133 Table<2,int> adjust_quad_dof_index_for_face_orientation_table;
02134
02156 std::vector<int> adjust_line_dof_index_for_line_orientation_table;
02157
02180 TableIndices<2>
02181 interface_constraints_size () const;
02182
02188 void compute_2nd (const Mapping<dim,spacedim> &mapping,
02189 const typename Triangulation<dim,spacedim>::cell_iterator &cell,
02190 const unsigned int offset,
02191 typename Mapping<dim,spacedim>::InternalDataBase &mapping_internal,
02192 InternalDataBase &fe_internal,
02193 FEValuesData<dim,spacedim> &data) const;
02194
02205 static
02206 std::vector<unsigned int>
02207 compute_n_nonzero_components (const std::vector<std::vector<bool> > &nonzero_components);
02208
02248 virtual UpdateFlags update_once (const UpdateFlags flags) const = 0;
02249
02265 virtual UpdateFlags update_each (const UpdateFlags flags) const = 0;
02266
02278 virtual FiniteElement<dim,spacedim> *clone() const = 0;
02279
02280 private:
02286 std::vector< std::pair<unsigned int, unsigned int> > system_to_component_table;
02287
02304 std::vector< std::pair<unsigned int, unsigned int> > face_system_to_component_table;
02305
02339 std::vector<std::pair<std::pair<unsigned int,unsigned int>,unsigned int> >
02340 system_to_base_table;
02341
02346 std::vector<std::pair<std::pair<unsigned int,unsigned int>,unsigned int> >
02347 face_system_to_base_table;
02348
02355 BlockIndices base_to_block_indices;
02356
02384 std::vector<std::pair<std::pair<unsigned int, unsigned int>, unsigned int> >
02385 component_to_base_table;
02386
02458 const std::vector<bool> restriction_is_additive_flags;
02459
02473 const std::vector<std::vector<bool> > nonzero_components;
02474
02489 const std::vector<unsigned int> n_nonzero_components_table;
02490
02502 static const double fd_step_length;
02503
02515 virtual typename Mapping<dim,spacedim>::InternalDataBase*
02516 get_data (const UpdateFlags flags,
02517 const Mapping<dim,spacedim> &mapping,
02518 const Quadrature<dim> &quadrature) const = 0;
02519
02532 virtual typename Mapping<dim,spacedim>::InternalDataBase*
02533 get_face_data (const UpdateFlags flags,
02534 const Mapping<dim,spacedim> &mapping,
02535 const Quadrature<dim-1> &quadrature) const;
02536
02549 virtual typename Mapping<dim,spacedim>::InternalDataBase*
02550 get_subface_data (const UpdateFlags flags,
02551 const Mapping<dim,spacedim> &mapping,
02552 const Quadrature<dim-1> &quadrature) const;
02553
02565 virtual void
02566 fill_fe_values (const Mapping<dim,spacedim> &mapping,
02567 const typename Triangulation<dim,spacedim>::cell_iterator &cell,
02568 const Quadrature<dim> &quadrature,
02569 typename Mapping<dim,spacedim>::InternalDataBase &mapping_internal,
02570 typename Mapping<dim,spacedim>::InternalDataBase &fe_internal,
02571 FEValuesData<dim,spacedim> &data,
02572 CellSimilarity::Similarity &cell_similarity) const = 0;
02573
02585 virtual void
02586 fill_fe_face_values (const Mapping<dim,spacedim> &mapping,
02587 const typename Triangulation<dim,spacedim>::cell_iterator &cell,
02588 const unsigned int face_no,
02589 const Quadrature<dim-1> &quadrature,
02590 typename Mapping<dim,spacedim>::InternalDataBase &mapping_internal,
02591 typename Mapping<dim,spacedim>::InternalDataBase &fe_internal,
02592 FEValuesData<dim,spacedim> &data) const = 0;
02593
02605 virtual void
02606 fill_fe_subface_values (const Mapping<dim,spacedim> &mapping,
02607 const typename Triangulation<dim,spacedim>::cell_iterator &cell,
02608 const unsigned int face_no,
02609 const unsigned int sub_no,
02610 const Quadrature<dim-1> &quadrature,
02611 typename Mapping<dim,spacedim>::InternalDataBase &mapping_internal,
02612 typename Mapping<dim,spacedim>::InternalDataBase &fe_internal,
02613 FEValuesData<dim,spacedim> &data) const = 0;
02614
02615 friend class InternalDataBase;
02616 friend class FEValuesBase<dim,spacedim>;
02617 friend class FEValues<dim,spacedim>;
02618 friend class FEFaceValues<dim,spacedim>;
02619 friend class FESubfaceValues<dim,spacedim>;
02620 template <int, int > friend class FESystem;
02621 template <class POLY, int dim_, int spacedim_> friend class FE_PolyTensor;
02622 friend class hp::FECollection<dim,spacedim>;
02623
02624 };
02625
02626
02627
02628
02629
02630 template <int dim, int spacedim>
02631 inline
02632 const FiniteElement<dim,spacedim> &
02633 FiniteElement<dim,spacedim>::operator[] (const unsigned int fe_index) const
02634 {
02635 Assert (fe_index == 0,
02636 ExcMessage ("A fe_index of zero is the only index allowed here"));
02637 return *this;
02638 }
02639
02640
02641
02642 template <int dim, int spacedim>
02643 inline
02644 std::pair<unsigned int,unsigned int>
02645 FiniteElement<dim,spacedim>::system_to_component_index (const unsigned int index) const
02646 {
02647 Assert (index < system_to_component_table.size(),
02648 ExcIndexRange(index, 0, system_to_component_table.size()));
02649 Assert (is_primitive (index),
02650 ( typename FiniteElement<dim,spacedim>::ExcShapeFunctionNotPrimitive(index)) );
02651 return system_to_component_table[index];
02652 }
02653
02654
02655
02656 template <int dim, int spacedim>
02657 inline
02658 unsigned int
02659 FiniteElement<dim,spacedim>::n_base_elements () const
02660 {
02661 return base_to_block_indices.size();
02662 }
02663
02664
02665
02666 template <int dim, int spacedim>
02667 inline
02668 unsigned int
02669 FiniteElement<dim,spacedim>::element_multiplicity (const unsigned int index) const
02670 {
02671 return base_to_block_indices.block_size(index);
02672 }
02673
02674
02675
02676 template <int dim, int spacedim>
02677 inline
02678 unsigned int
02679 FiniteElement<dim,spacedim>::component_to_system_index (const unsigned int component,
02680 const unsigned int index) const
02681 {
02682 std::vector< std::pair<unsigned int, unsigned int> >::const_iterator
02683 it = std::find(system_to_component_table.begin(), system_to_component_table.end(),
02684 std::pair<unsigned int, unsigned int>(component, index));
02685
02686 Assert(it != system_to_component_table.end(), ExcComponentIndexInvalid(component, index));
02687 return std::distance(system_to_component_table.begin(), it);
02688 }
02689
02690
02691
02692 template <int dim, int spacedim>
02693 inline
02694 std::pair<unsigned int,unsigned int>
02695 FiniteElement<dim,spacedim>::face_system_to_component_index (const unsigned int index) const
02696 {
02697 Assert(index < face_system_to_component_table.size(),
02698 ExcIndexRange(index, 0, face_system_to_component_table.size()));
02699
02700
02701
02702
02703
02704
02705
02706
02707
02708
02709
02710
02711
02712
02713 Assert (is_primitive(this->face_to_equivalent_cell_index(index)),
02714 (typename FiniteElement<dim,spacedim>::ExcShapeFunctionNotPrimitive(index)) );
02715
02716 return face_system_to_component_table[index];
02717 }
02718
02719
02720
02721 template <int dim, int spacedim>
02722 inline
02723 std::pair<std::pair<unsigned int,unsigned int>,unsigned int>
02724 FiniteElement<dim,spacedim>::system_to_base_index (const unsigned int index) const
02725 {
02726 Assert (index < system_to_base_table.size(),
02727 ExcIndexRange(index, 0, system_to_base_table.size()));
02728 return system_to_base_table[index];
02729 }
02730
02731
02732
02733
02734 template <int dim, int spacedim>
02735 inline
02736 std::pair<std::pair<unsigned int,unsigned int>,unsigned int>
02737 FiniteElement<dim,spacedim>::face_system_to_base_index (const unsigned int index) const
02738 {
02739 Assert(index < face_system_to_base_table.size(),
02740 ExcIndexRange(index, 0, face_system_to_base_table.size()));
02741 return face_system_to_base_table[index];
02742 }
02743
02744
02745
02746 template <int dim, int spacedim>
02747 inline
02748 unsigned int
02749 FiniteElement<dim,spacedim>::first_block_of_base (const unsigned int index) const
02750 {
02751 return base_to_block_indices.block_start(index);
02752 }
02753
02754
02755
02756 template <int dim, int spacedim>
02757 inline
02758 std::pair<unsigned int,unsigned int>
02759 FiniteElement<dim,spacedim>::component_to_base_index (const unsigned int index) const
02760 {
02761 Assert(index < component_to_base_table.size(),
02762 ExcIndexRange(index, 0, component_to_base_table.size()));
02763
02764 return component_to_base_table[index].first;
02765 }
02766
02767
02768
02769 template <int dim, int spacedim>
02770 inline
02771 std::pair<unsigned int,unsigned int>
02772 FiniteElement<dim,spacedim>::block_to_base_index (const unsigned int index) const
02773 {
02774 return base_to_block_indices.global_to_local(index);
02775 }
02776
02777
02778
02779 template <int dim, int spacedim>
02780 inline
02781 std::pair<unsigned int,unsigned int>
02782 FiniteElement<dim,spacedim>::system_to_block_index (const unsigned int index) const
02783 {
02784 Assert (index < this->dofs_per_cell,
02785 ExcIndexRange(index, 0, this->dofs_per_cell));
02786
02787
02788
02789 return std::pair<unsigned int, unsigned int>(
02790 first_block_of_base(system_to_base_table[index].first.first)
02791 + system_to_base_table[index].first.second,
02792 system_to_base_table[index].second);
02793 }
02794
02795
02796
02797 template <int dim, int spacedim>
02798 inline
02799 bool
02800 FiniteElement<dim,spacedim>::restriction_is_additive (const unsigned int index) const
02801 {
02802 Assert(index < this->dofs_per_cell,
02803 ExcIndexRange(index, 0, this->dofs_per_cell));
02804 return restriction_is_additive_flags[index];
02805 }
02806
02807
02808
02809 template <int dim, int spacedim>
02810 inline
02811 const std::vector<bool> &
02812 FiniteElement<dim,spacedim>::get_nonzero_components (const unsigned int i) const
02813 {
02814 Assert (i < this->dofs_per_cell, ExcIndexRange (i, 0, this->dofs_per_cell));
02815 return nonzero_components[i];
02816 }
02817
02818
02819
02820 template <int dim, int spacedim>
02821 inline
02822 unsigned int
02823 FiniteElement<dim,spacedim>::n_nonzero_components (const unsigned int i) const
02824 {
02825 Assert (i < this->dofs_per_cell, ExcIndexRange (i, 0, this->dofs_per_cell));
02826 return n_nonzero_components_table[i];
02827 }
02828
02829
02830
02831 template <int dim, int spacedim>
02832 inline
02833 bool
02834 FiniteElement<dim,spacedim>::is_primitive (const unsigned int i) const
02835 {
02836 Assert (i < this->dofs_per_cell, ExcIndexRange (i, 0, this->dofs_per_cell));
02837
02838
02839
02840
02841
02842
02843
02844
02845
02846
02847
02848
02849
02850 return (is_primitive() ||
02851 (n_nonzero_components_table[i] == 1));
02852 }
02853
02854
02855
02856 DEAL_II_NAMESPACE_CLOSE
02857
02858 #endif
02859