mymapping_q1.h

Go to the documentation of this file.

#ifndef _deal2__mymapping_q1_h
#define _deal2__mymapping_q1_h


#include <base/config.h>
#include <base/table.h>
#include <base/qprojector.h>
#include <grid/tria_iterator.h>
#include <dofs/dof_accessor.h>
#include <fe/mapping.h>

#include <cmath>
#include <iostream>

using namespace dealii;

template <int dim>
class MyMappingQ1 : public Mapping<dim>
{
  public:
    MyMappingQ1 ();

    virtual Point<dim>
      transform_unit_to_real_cell (
          const typename Triangulation<dim>::cell_iterator &cell,
          const Point<dim>                                 &p) const;

    virtual Point<dim>
      transform_real_to_unit_cell (
          const typename Triangulation<dim>::cell_iterator &cell,
          const Point<dim>                                 &p) const;
          

    virtual void
      transform_covariant (const VectorSlice<const std::vector<Tensor<1,dim> > > input,
          const unsigned int                 offset,
          VectorSlice<std::vector<Tensor<1,dim> > > output,
          const typename Mapping<dim>::InternalDataBase &internal) const;
          
     
    virtual void
      transform_covariant (const VectorSlice<const std::vector<Tensor<2,dim> > > input,
          const unsigned int                 offset,
          VectorSlice<std::vector<Tensor<2,dim> > > output,
          const typename Mapping<dim>::InternalDataBase &internal) const;
          
    
    virtual void
      transform_contravariant (const VectorSlice<const std::vector<Tensor<1,dim> > > input,
          const unsigned int                 offset,
          VectorSlice<std::vector<Tensor<1,dim> > > output,
          const typename Mapping<dim>::InternalDataBase &internal) const;
          

    virtual void
      transform_contravariant (const VectorSlice<const std::vector<Tensor<2,dim> > > input,
          const unsigned int                 offset,
          VectorSlice<std::vector<Tensor<2,dim> > > output,
          const typename Mapping<dim>::InternalDataBase &internal) const;
          
    
    virtual UpdateFlags update_once (const UpdateFlags flags) const;
     

    virtual UpdateFlags update_each (const UpdateFlags flags) const;

    virtual
      Mapping<dim> * clone () const;
      
    class InternalData : public Mapping<dim>::InternalDataBase
  {
    public:
      InternalData(const unsigned int n_shape_functions);

      double shape (const unsigned int qpoint,
          const unsigned int shape_nr) const;

      double &shape (const unsigned int qpoint,
          const unsigned int shape_nr);

      double support_shape (const unsigned int point_no,
          const unsigned int shape_nr) const;

      double &support_shape (const unsigned int point_no,
          const unsigned int shape_nr);

      Tensor<1,dim> derivative (const unsigned int qpoint,
          const unsigned int shape_nr) const;

      Tensor<1,dim> &derivative (const unsigned int qpoint,
          const unsigned int shape_nr);

      Tensor<1,dim> support_derivative (const unsigned int point_no,
          const unsigned int shape_nr) const;

      Tensor<1,dim> &support_derivative (const unsigned int point_no,
          const unsigned int shape_nr);

      virtual unsigned int memory_consumption () const;

      std::vector<double> shape_values;

      std::vector<Tensor<1,dim> > shape_derivatives;

      std::vector<Tensor<2,dim> > covariant;

      std::vector<Tensor<2,dim> > contravariant;

      std::vector<double> support_values;

      std::vector<Tensor<1,dim> > support_derivatives;

      std::vector<std::vector<Tensor<1,dim> > > unit_tangentials;

      std::vector<std::vector<Tensor<1,dim> > > aux;

      std::vector<Point<dim> > mapping_support_points;

      std::vector<Point<dim> > mapping_generalized_support_points;

      typename Triangulation<dim>::cell_iterator cell_of_current_support_points;

      bool is_mapping_q1_data;

      unsigned int n_shape_functions;
  };

    DeclException0 (ExcAccessToUninitializedField);

  protected:

    typedef
      typename QProjector<dim>::DataSetDescriptor
      DataSetDescriptor;

  //virtual void
  //  fill_support_points (
  //      const typename Triangulation<dim>::cell_iterator &cell,
  //      typename Mapping<dim>::InternalDataBase          &mapping_data,
  //      const std::vector<Point<dim> >                   &generalized_support_points) const ;


  virtual void
    fill_fe_values (
        const typename Triangulation<dim>::cell_iterator &cell,
        const Quadrature<dim>                            &quadrature,
        typename Mapping<dim>::InternalDataBase          &mapping_data,
        typename std::vector<Point<dim> >                &quadrature_points,
        std::vector<double>                              &JxW_values,
        std::vector<Tensor<2,dim> >                      &jacobians, 
        std::vector<Tensor<3,dim> >                      &jacobian_grads,
        std::vector< Tensor< 2, dim > >                  &inverse_jacobians) const;

    virtual void
      fill_fe_face_values (const typename Triangulation<dim>::cell_iterator &cell,
          const unsigned int face_no,
          const Quadrature<dim-1>& quadrature,
          typename Mapping<dim>::InternalDataBase &mapping_data,
          typename std::vector<Point<dim> >        &quadrature_points,
          std::vector<double>             &JxW_values,
          typename std::vector<Tensor<1,dim> >        &boundary_form,
          typename std::vector<Point<dim> >        &normal_vectors,
          std::vector<double>             &cell_JxW_values) const ;

    virtual void
      fill_fe_subface_values (const typename Triangulation<dim>::cell_iterator &cell,
          const unsigned int face_no,
          const unsigned int sub_no,
          const Quadrature<dim-1>& quadrature,
          typename Mapping<dim>::InternalDataBase &mapping_data,
          typename std::vector<Point<dim> >        &quadrature_points,
          std::vector<double>             &JxW_values,
          typename std::vector<Tensor<1,dim> >        &boundary_form,
          typename std::vector<Point<dim> >        &normal_vectors,
          std::vector<double>             &cell_JxW_values) const ;

    void compute_shapes (
        const std::vector<Point<dim> > &quadrature_points,
        InternalData &data) const;

    void compute_data (
        const UpdateFlags               flags,
        const std::vector<Point<dim> > &quadrature_points,
        const unsigned int              n_orig_q_points,
        InternalData                   &data) const;

    //void compute_data (
    //    const UpdateFlags               flags,
    //    const std::vector<Point<dim> > &generalized_support_points,
    //    const std::vector<Point<dim> > &quadrature_points,
    //    const unsigned int              n_orig_q_points,
    //    InternalData                   &data) const;

    void compute_face_data (
        const UpdateFlags flags,
        const Quadrature<dim> &quadrature,
        const unsigned int n_orig_q_points,
        InternalData &data) const;

    void compute_fill (const typename Triangulation<dim>::cell_iterator &cell,
        const unsigned int      npts,
        const DataSetDescriptor data_set,
        InternalData           &data,
        std::vector<Point<dim> > &quadrature_points) const;

    void compute_fill_face (const typename Triangulation<dim>::cell_iterator &cell,
        const unsigned int      face_no,
        const bool              is_subface,
        const unsigned int      npts,
        const DataSetDescriptor data_set,
        const std::vector<double>   &weights,
        InternalData           &mapping_data,
        std::vector<Point<dim> >    &quadrature_points,
        std::vector<double>         &JxW_values,
        std::vector<Tensor<1,dim> > &boundary_form,
        std::vector<Point<dim> > &normal_vectors,
        std::vector<double>         &cell_JxW_values) const;

    virtual void compute_shapes_virtual (
        const std::vector<Point<dim> > &quadrature_points,
        InternalData &data) const;

    Point<dim> transform_unit_to_real_cell_internal (const InternalData &mdata) const;

    void transform_real_to_unit_cell_internal (const typename Triangulation<dim>::cell_iterator &cell,
        const Point<dim> &p,
        InternalData &mdata,
        Point<dim> &p_unit) const;

  private:
    virtual
      typename Mapping<dim>::InternalDataBase *
      get_data (
          const UpdateFlags,
          const Quadrature<dim>          & quadrature) const;

    //virtual
    //  typename Mapping<dim>::InternalDataBase *
    //  get_data (
    //      const UpdateFlags,
    //      const Quadrature<dim>          & quadrature,
    //      const std::vector<Point<dim> > & points) const;

    virtual
      typename Mapping<dim>::InternalDataBase *
      get_face_data (const UpdateFlags flags,
          const Quadrature<dim-1>& quadrature) const;

    virtual
      typename Mapping<dim>::InternalDataBase *
      get_subface_data (const UpdateFlags flags,
          const Quadrature<dim-1>& quadrature) const;

    virtual void compute_mapping_support_points(
        const typename Triangulation<dim>::cell_iterator &cell,
        std::vector<Point<dim> > &a) const;

    static const unsigned int n_shape_functions = GeometryInfo<dim>::vertices_per_cell;
};


//template <int dim>
//struct StaticMappingQ1
//{
//  static MyMappingQ1<dim> mapping;
//};



/*----------------------------------------------------------------------*/

template<int dim>
inline
double
MyMappingQ1<dim>::InternalData::shape (const unsigned int qpoint,
    const unsigned int shape_nr) const
{
  Assert(qpoint*n_shape_functions + shape_nr < shape_values.size(),
      ExcIndexRange(qpoint*n_shape_functions + shape_nr, 0,
        shape_values.size()));
  return shape_values [qpoint*n_shape_functions + shape_nr];
}



template <int dim>
inline
  double &
MyMappingQ1<dim>::InternalData::shape (const unsigned int qpoint,
    const unsigned int shape_nr)
{
  Assert(qpoint*n_shape_functions + shape_nr < shape_values.size(),
      ExcIndexRange(qpoint*n_shape_functions + shape_nr, 0,
        shape_values.size()));
  return shape_values [qpoint*n_shape_functions + shape_nr];
}

template<int dim>
inline
double
MyMappingQ1<dim>::InternalData::support_shape (const unsigned int point_no,
    const unsigned int shape_nr) const
{
  Assert(point_no*n_shape_functions + shape_nr < support_values.size(),
      ExcIndexRange(point_no*n_shape_functions + shape_nr, 0,
        support_values.size()));
  return support_values [point_no*n_shape_functions + shape_nr];
}



template <int dim>
inline
  double &
MyMappingQ1<dim>::InternalData::support_shape (const unsigned int point_no,
    const unsigned int shape_nr)
{
  Assert(point_no*n_shape_functions + shape_nr < support_values.size(),
      ExcIndexRange(point_no*n_shape_functions + shape_nr, 0,
        support_values.size()));
  return support_values [point_no*n_shape_functions + shape_nr];
}

template <int dim>
inline
Tensor<1,dim>
MyMappingQ1<dim>::InternalData::derivative (const unsigned int qpoint,
    const unsigned int shape_nr) const
{
  Assert(qpoint*n_shape_functions + shape_nr < shape_derivatives.size(),
      ExcIndexRange(qpoint*n_shape_functions + shape_nr, 0,
        shape_derivatives.size()));
  return shape_derivatives [qpoint*n_shape_functions + shape_nr];
}

template <int dim>
inline
  Tensor<1,dim> &
MyMappingQ1<dim>::InternalData::derivative (const unsigned int qpoint,
    const unsigned int shape_nr)
{
  Assert(qpoint*n_shape_functions + shape_nr < shape_derivatives.size(),
      ExcIndexRange(qpoint*n_shape_functions + shape_nr, 0,
        shape_derivatives.size()));
  return shape_derivatives [qpoint*n_shape_functions + shape_nr];
}

template <int dim>
inline
Tensor<1,dim>
MyMappingQ1<dim>::InternalData::support_derivative (const unsigned int point_no,
    const unsigned int shape_nr) const
{
  Assert(point_no*n_shape_functions + shape_nr < support_derivatives.size(),
      ExcIndexRange(point_no*n_shape_functions + shape_nr, 0,
        support_derivatives.size()));
  return support_derivatives [point_no*n_shape_functions + shape_nr];
}

template <int dim>
inline
  Tensor<1,dim> &
MyMappingQ1<dim>::InternalData::support_derivative (const unsigned int point_no,
    const unsigned int shape_nr)
{
  Assert(point_no*n_shape_functions + shape_nr < support_derivatives.size(),
      ExcIndexRange(point_no*n_shape_functions + shape_nr, 0,
        support_derivatives.size()));
  return support_derivatives [point_no*n_shape_functions + shape_nr];
}


/* -------------- declaration of explicit specializations ------------- */

#ifndef DOXYGEN

template<> void MyMappingQ1<1>::compute_shapes_virtual (
    const std::vector<Point<1> > &quadrature_points,
    InternalData& data) const;
template<> void MyMappingQ1<2>::compute_shapes_virtual (
    const std::vector<Point<2> > &quadrature_points,
    InternalData &data) const;
template<> void MyMappingQ1<3>::compute_shapes_virtual (
    const std::vector<Point<3> > &quadrature_points,
    InternalData &data) const;

template <>
void
MyMappingQ1<1>::compute_face_data (const UpdateFlags,
    const Quadrature<1> &,
    const unsigned int,
    InternalData &) const;

template <> void MyMappingQ1<1>::compute_fill_face (
    const Triangulation<1>::cell_iterator &,
    const unsigned int,
    const bool,
    const unsigned int,
    const DataSetDescriptor,
    const std::vector<double> &,
    InternalData &,
    std::vector<Point<1> > &,
    std::vector<double> &,
    std::vector<Tensor<1,1> > &,
    std::vector<Point<1> > &,
    std::vector<double> &) const;

template <> void MyMappingQ1<1>::fill_fe_face_values (
    const Triangulation<1>::cell_iterator &,
    const unsigned,
    const Quadrature<0>&,
    Mapping<1>::InternalDataBase&,
    std::vector<Point<1> >&,
    std::vector<double>&,
    std::vector<Tensor<1,1> >&,
    std::vector<Point<1> >&,
    std::vector<double>&) const;

template <> void MyMappingQ1<1>::fill_fe_subface_values (
    const Triangulation<1>::cell_iterator &,
    const unsigned,
    const unsigned,
    const Quadrature<0>&,
    Mapping<1>::InternalDataBase&,
    std::vector<Point<1> >&,
    std::vector<double>&,
    std::vector<Tensor<1,1> >&,
    std::vector<Point<1> >&,
    std::vector<double>&) const;

#endif // DOXYGEN

#endif // _deal2__mymapping_q1_h

---