fem_context.h

Go to the documentation of this file.

// The libMesh Finite Element Library.
// Copyright (C) 2002-2014 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner

// This library is free software; you can redistribute it and/or
// modify it under the terms of the GNU Lesser General Public
// License as published by the Free Software Foundation; either
// version 2.1 of the License, or (at your option) any later version.

// This library is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
// Lesser General Public License for more details.

// You should have received a copy of the GNU Lesser General Public
// License along with this library; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA



#ifndef LIBMESH_FEM_CONTEXT_H
#define LIBMESH_FEM_CONTEXT_H

// Local Includes
#include "libmesh/diff_context.h"
#include "libmesh/id_types.h"
#include "libmesh/fe_type.h"
#include "libmesh/fe_base.h"
#include "libmesh/vector_value.h"

#ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
#include "libmesh/tensor_value.h"
#endif

// C++ includes
#include <map>

namespace libMesh
{

// Forward Declarations
class BoundaryInfo;
class Elem;
template <typename T> class FEGenericBase;
typedef FEGenericBase<Real> FEBase;
class QBase;
class Point;
template <typename T> class NumericVector;

// ------------------------------------------------------------
// FEMContext class definition

class FEMContext : public DiffContext
{
public:

  explicit
  FEMContext (const System &sys);

  virtual ~FEMContext ();

  bool has_side_boundary_id(boundary_id_type id) const;

  std::vector<boundary_id_type> side_boundary_ids() const;

  Number interior_value(unsigned int var, unsigned int qp) const;

  Number side_value(unsigned int var, unsigned int qp) const;

  Number point_value(unsigned int var, const Point &p) const;

  Gradient interior_gradient(unsigned int var, unsigned int qp) const;

  Gradient side_gradient(unsigned int var, unsigned int qp) const;

  Gradient point_gradient(unsigned int var, const Point &p) const;

#ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES

  Tensor interior_hessian(unsigned int var, unsigned int qp) const;

  Tensor side_hessian(unsigned int var, unsigned int qp) const;

  Tensor point_hessian(unsigned int var, const Point &p) const;

#endif // LIBMESH_ENABLE_SECOND_DERIVATIVES

  Number fixed_interior_value(unsigned int var, unsigned int qp) const;

  Number fixed_side_value(unsigned int var, unsigned int qp) const;

  Number fixed_point_value(unsigned int var, const Point &p) const;

  Gradient fixed_interior_gradient(unsigned int var, unsigned int qp) const;

  Gradient fixed_side_gradient(unsigned int var, unsigned int qp) const;

  Gradient fixed_point_gradient(unsigned int var, const Point &p) const;

#ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES

  Tensor fixed_interior_hessian(unsigned int var, unsigned int qp) const;

  Tensor fixed_side_hessian(unsigned int var, unsigned int qp) const;

  Tensor fixed_point_hessian (unsigned int var, const Point &p) const;

#endif // LIBMESH_ENABLE_SECOND_DERIVATIVES

  template<typename OutputShape>
  void get_element_fe( unsigned int var, FEGenericBase<OutputShape> *& fe ) const
  { this->get_element_fe<OutputShape>(var,fe,this->get_dim()); }

  FEBase* get_element_fe( unsigned int var ) const
  { return this->get_element_fe(var,this->get_dim()); }

  template<typename OutputShape>
  void get_element_fe( unsigned int var, FEGenericBase<OutputShape> *& fe,
                       unsigned char dim ) const;

  FEBase* get_element_fe( unsigned int var, unsigned char dim ) const;

  template<typename OutputShape>
  void get_side_fe( unsigned int var, FEGenericBase<OutputShape> *& fe ) const
  { this->get_side_fe<OutputShape>(var,fe,this->get_dim()); }

  FEBase* get_side_fe( unsigned int var ) const
  { return this->get_side_fe(var,this->get_dim()); }

  template<typename OutputShape>
  void get_side_fe( unsigned int var, FEGenericBase<OutputShape> *& fe,
                    unsigned char dim ) const;

  FEBase* get_side_fe( unsigned int var, unsigned char dim ) const;

  template<typename OutputShape>
  void get_edge_fe( unsigned int var, FEGenericBase<OutputShape> *& fe ) const;

  FEBase* get_edge_fe( unsigned int var ) const;

  template<typename OutputType>
  void interior_value(unsigned int var, unsigned int qp,
                      OutputType& u) const;

  template<typename OutputType>
  void interior_values(unsigned int var, const NumericVector<Number> & _system_vector,
                       std::vector<OutputType>& interior_values_vector) const;

  template<typename OutputType>
  void side_value(unsigned int var, unsigned int qp,
                  OutputType& u) const;

  template<typename OutputType>
  void side_values(unsigned int var, const NumericVector<Number> & _system_vector,
                   std::vector<OutputType>& side_values_vector) const;

  template<typename OutputType>
  void point_value(unsigned int var, const Point &p,
                   OutputType& u) const;

  template<typename OutputType>
  void interior_gradient(unsigned int var, unsigned int qp,
                         OutputType& du) const;

  template<typename OutputType>
  void interior_gradients(unsigned int var, const NumericVector<Number> & _system_vector,
                          std::vector<OutputType>& interior_gradients_vector) const;

  template<typename OutputType>
  void side_gradient(unsigned int var, unsigned int qp,
                     OutputType & du) const;

  template<typename OutputType>
  void side_gradients(unsigned int var, const NumericVector<Number> & _system_vector,
                      std::vector<OutputType>& side_gradients_vector) const;

  template<typename OutputType>
  void point_gradient(unsigned int var, const Point &p,
                      OutputType& grad_u) const;

#ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES

  template<typename OutputType>
  void interior_hessian(unsigned int var, unsigned int qp,
                        OutputType& d2u) const;

  template<typename OutputType>
  void interior_hessians(unsigned int var, const NumericVector<Number> & _system_vector,
                         std::vector<OutputType>& d2u_vals) const;

  template<typename OutputType>
  void side_hessian(unsigned int var, unsigned int qp,
                    OutputType& d2u) const;

  template<typename OutputType>
  void side_hessians(unsigned int var, const NumericVector<Number> & _system_vector,
                     std::vector<OutputType>& d2u_vals) const;

  template<typename OutputType>
  void point_hessian(unsigned int var, const Point &p,
                     OutputType& hess_u) const;

#endif // LIBMESH_ENABLE_SECOND_DERIVATIVES

  template<typename OutputType>
  void interior_rate(unsigned int var, unsigned int qp,
                     OutputType& u) const;

  template<typename OutputType>
  void side_rate(unsigned int var, unsigned int qp,
                 OutputType& u) const;

  template<typename OutputType>
  void point_rate(unsigned int var, const Point &p,
                  OutputType& u) const;

  template<typename OutputType>
  void fixed_interior_value(unsigned int var, unsigned int qp,
                            OutputType& u) const;

  template<typename OutputType>
  void fixed_side_value(unsigned int var, unsigned int qp,
                        OutputType& u) const;

  template<typename OutputType>
  void fixed_point_value(unsigned int var, const Point &p,
                         OutputType& u) const;

  template<typename OutputType>
  void fixed_interior_gradient(unsigned int var, unsigned int qp,
                               OutputType& grad_u) const;

  template<typename OutputType>
  void fixed_side_gradient(unsigned int var, unsigned int qp,
                           OutputType& grad_u) const;

  template<typename OutputType>
  void fixed_point_gradient(unsigned int var, const Point &p,
                            OutputType& grad_u) const;

#ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES

  template<typename OutputType>
  void fixed_interior_hessian(unsigned int var, unsigned int qp,
                              OutputType& hess_u) const;

  template<typename OutputType>
  void fixed_side_hessian(unsigned int var, unsigned int qp,
                          OutputType& hess_u) const;

  template<typename OutputType>
  void fixed_point_hessian(unsigned int var, const Point &p,
                           OutputType& hess_u) const;

#endif // LIBMESH_ENABLE_SECOND_DERIVATIVES

  template<typename OutputType>
  void interior_curl(unsigned int var, unsigned int qp,
                     OutputType& curl_u) const;

  template<typename OutputType>
  void point_curl(unsigned int var, const Point &p,
                  OutputType& curl_u) const;

  template<typename OutputType>
  void interior_div(unsigned int var, unsigned int qp,
                    OutputType& div_u) const;

  //  virtual void reinit(const FEMSystem&, Elem*);

  // should be protected:
  virtual void elem_reinit(Real theta);

  virtual void elem_side_reinit(Real theta);

  virtual void elem_edge_reinit(Real theta);

  virtual void nonlocal_reinit(Real theta);

  virtual void pre_fe_reinit(const System&, const Elem *e);

  void elem_fe_reinit();

  virtual void side_fe_reinit();

  void edge_fe_reinit();

  const QBase& get_element_qrule() const
  { return this->get_element_qrule(this->get_elem_dim()); }

  const QBase& get_side_qrule() const
  { return this->get_side_qrule(this->get_elem_dim()); }

  const QBase& get_element_qrule( unsigned char dim ) const
  { libmesh_assert(_element_qrule[dim]);
    return *(this->_element_qrule[dim]); }

  const QBase& get_side_qrule( unsigned char dim ) const
  { libmesh_assert(_side_qrule[dim]);
    return *(this->_side_qrule[dim]); }

  const QBase& get_edge_qrule() const
  { return *(this->_edge_qrule); }

  virtual void set_mesh_system(System* sys)
  { this->_mesh_sys = sys; }

  const System* get_mesh_system() const
  { return this->_mesh_sys; }

  System* get_mesh_system()
  { return this->_mesh_sys; }

  unsigned int get_mesh_x_var() const
  { return _mesh_x_var; }

  void set_mesh_x_var(unsigned int x_var)
  { _mesh_x_var = x_var; }

  unsigned int get_mesh_y_var() const
  { return _mesh_y_var; }

  void set_mesh_y_var(unsigned int y_var)
  { _mesh_y_var = y_var; }

  unsigned int get_mesh_z_var() const
  { return _mesh_z_var; }

  void set_mesh_z_var(unsigned int z_var)
  { _mesh_z_var = z_var; }

  bool has_elem() const
  { return (this->_elem != NULL); }

  const Elem& get_elem() const
  { libmesh_assert(this->_elem);
    return *(this->_elem); }

  Elem& get_elem()
  { libmesh_assert(this->_elem);
    return *(const_cast<Elem*>(this->_elem)); }

  unsigned char get_side() const
  { return side; }

  unsigned char get_edge() const
  { return edge; }

  unsigned char get_dim() const
  { return this->_dim; }

  unsigned char get_elem_dim() const
  { return _elem_dim; }


  void elem_position_set(Real theta);

  void elem_position_get();

  System *_mesh_sys;

  unsigned int _mesh_x_var, _mesh_y_var, _mesh_z_var;

  unsigned char side;

  unsigned char edge;

protected:

  template<typename OutputShape>
  UniquePtr<FEGenericBase<OutputShape> > build_new_fe( const FEGenericBase<OutputShape>* fe, const Point &p ) const;

  void set_elem( const Elem* e );

  // gcc-3.4, oracle 12.3 require this typedef to be public
  // in order to use it in a return type
public:

  typedef const DenseSubVector<Number>& (DiffContext::*diff_subsolution_getter)(unsigned int) const;

protected:
  template<typename OutputType,
           diff_subsolution_getter subsolution_getter>
  void some_interior_value(unsigned int var, unsigned int qp, OutputType& u) const;

  template<typename OutputType,
           diff_subsolution_getter subsolution_getter>
  void some_interior_gradient(unsigned int var, unsigned int qp, OutputType& u) const;

#ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES

  template<typename OutputType,
           diff_subsolution_getter subsolution_getter>
  void some_interior_hessian(unsigned int var, unsigned int qp, OutputType& u) const;
#endif

  template<typename OutputType,
           diff_subsolution_getter subsolution_getter>
  void some_side_value(unsigned int var, unsigned int qp, OutputType& u) const;


  std::vector<std::map<FEType, FEAbstract*> > _element_fe;
  std::vector<std::map<FEType, FEAbstract*> > _side_fe;
  std::map<FEType, FEAbstract*> _edge_fe;


  std::vector<std::vector<FEAbstract*> > _element_fe_var;
  std::vector<std::vector<FEAbstract*> > _side_fe_var;
  std::vector<FEAbstract*> _edge_fe_var;

  const BoundaryInfo& _boundary_info;

  const Elem *_elem;

  unsigned char _dim;

  unsigned char _elem_dim;

  std::vector<QBase*> _element_qrule;

  std::vector<QBase*> _side_qrule;

  QBase *_edge_qrule;

private:
  void _do_elem_position_set(Real theta);

  void _update_time_from_system(Real theta);
};



// ------------------------------------------------------------
// FEMContext inline methods

inline
void FEMContext::elem_position_set(Real theta)
{
  if (_mesh_sys)
    this->_do_elem_position_set(theta);
}

template<typename OutputShape>
inline
void FEMContext::get_element_fe( unsigned int var, FEGenericBase<OutputShape> *& fe,
                                 unsigned char dim ) const
{
  libmesh_assert( !_element_fe_var[dim].empty() );
  libmesh_assert_less ( var, (_element_fe_var[dim].size() ) );
  fe = cast_ptr<FEGenericBase<OutputShape>*>( (_element_fe_var[dim][var] ) );
}

inline
FEBase* FEMContext::get_element_fe( unsigned int var, unsigned char dim ) const
{
  libmesh_assert( !_element_fe_var[dim].empty() );
  libmesh_assert_less ( var, (_element_fe_var[dim].size() ) );
  return cast_ptr<FEBase*>( (_element_fe_var[dim][var] ) );
}

template<typename OutputShape>
inline
void FEMContext::get_side_fe( unsigned int var, FEGenericBase<OutputShape> *& fe,
                              unsigned char dim ) const
{
  libmesh_assert( !_side_fe_var[dim].empty() );
  libmesh_assert_less ( var, (_side_fe_var[dim].size() ) );
  fe = cast_ptr<FEGenericBase<OutputShape>*>( (_side_fe_var[dim][var] ) );
}

inline
FEBase* FEMContext::get_side_fe( unsigned int var, unsigned char dim ) const
{
  libmesh_assert( !_side_fe_var[dim].empty() );
  libmesh_assert_less ( var, (_side_fe_var[dim].size() ) );
  return cast_ptr<FEBase*>( (_side_fe_var[dim][var] ) );
}

template<typename OutputShape>
inline
void FEMContext::get_edge_fe( unsigned int var, FEGenericBase<OutputShape> *& fe ) const
{
  libmesh_assert_less ( var, _edge_fe_var.size() );
  fe = cast_ptr<FEGenericBase<OutputShape>*>( _edge_fe_var[var] );
}

inline
FEBase* FEMContext::get_edge_fe( unsigned int var ) const
{
  libmesh_assert_less ( var, _edge_fe_var.size() );
  return cast_ptr<FEBase*>( _edge_fe_var[var] );
}


} // namespace libMesh

#endif // LIBMESH_FEM_CONTEXT_H