include/deal.II/fe/fe.h

00001 //---------------------------------------------------------------------------
00002 //    @f$Id: fe.h 25345 2012-03-31 08:37:04Z bangerth @f$
00003 //
00004 //    Copyright (C) 1998, 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2010, 2011, 2012 by the deal.II authors
00005 //
00006 //    This file is subject to QPL and may not be  distributed
00007 //    without copyright and license information. Please refer
00008 //    to the file deal.II/doc/license.html for the  text  and
00009 //    further information on this license.
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                                    // in debug mode, check whether the
02701                                    // function is primitive, since
02702                                    // otherwise the result may have no
02703                                    // meaning
02704                                    //
02705                                    // since the primitivity tables are
02706                                    // all geared towards cell dof
02707                                    // indices, rather than face dof
02708                                    // indices, we have to work a
02709                                    // little bit...
02710                                    //
02711                                    // in 1d, the face index is equal
02712                                    // to the cell index
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                                    // The block is computed simply as
02787                                    // first block of this base plus
02788                                    // the index within the base blocks
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                                    // return primitivity of a shape
02839                                    // function by checking whether it
02840                                    // has more than one non-zero
02841                                    // component or not. we could cache
02842                                    // this value in an array of bools,
02843                                    // but accessing a bit-vector (as
02844                                    // std::vector<bool> is) is
02845                                    // probably more expensive than
02846                                    // just comparing against 1
02847                                    //
02848                                    // for good measure, short circuit the test
02849                                    // if the entire FE is primitive
02850   return (is_primitive() ||
02851           (n_nonzero_components_table[i] == 1));
02852 }
02853 
02854 
02855 
02856 DEAL_II_NAMESPACE_CLOSE
02857 
02858 #endif
02859 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines

deal.II documentation generated on Tue May 22 2012 12:06:02 by doxygen 1.7.3