old_interface.h

Go to the documentation of this file.

#ifndef _interface_h_
#define _interface_h_

#include <deal.II/lac/linear_operator.h>

#include <deal2lkit/dof_utilities.h>
#include <deal2lkit/parsed_finite_element.h>
#include <deal2lkit/any_data.h>
#include <deal2lkit/parsed_function.h>
#include <deal2lkit/parsed_mapped_functions.h>
#include <deal2lkit/parsed_dirichlet_bcs.h>

#include "data/assembly.h"
#include "lac/lac_type.h"
#include <deal.II/base/sacado_product_type.h>

template <int dim,int spacedim=dim, int n_components=1, typename LAC=LATrilinos>
class Interface : public ParsedFiniteElement<dim,spacedim>
{
  typedef FEValuesCache<dim,spacedim> Scratch;
  typedef Assembly::CopyData::piDoMUSPreconditioner<dim,dim> CopyPreconditioner;
  typedef Assembly::CopyData::piDoMUSSystem<dim,dim> CopySystem;

public:

  virtual ~Interface() {};

  Interface(const std::string &name="",
            const std::string &default_fe="FE_Q(1)",
            const std::string &default_component_names="u",
            const std::string &default_coupling="",
            const std::string &default_preconditioner_coupling="",
            const std::string &default_differential_components="") :
    ParsedFiniteElement<dim,spacedim>(name, default_fe, default_component_names,
                                      n_components, default_coupling, default_preconditioner_coupling),
    forcing_terms("Forcing terms", default_component_names, "0=ALL"),
    neumann_bcs("Neumann boundary conditions", default_component_names, "0=ALL"),
    dirichlet_bcs("Dirichlet boundary conditions", default_component_names, "0=ALL"),
    str_diff_comp(default_differential_components),
    old_t(-1.0)
  {};

  virtual void declare_parameters (ParameterHandler &prm);
  virtual void parse_parameters_call_back ();
  const std::vector<unsigned int> get_differential_blocks() const;

  virtual void set_time (const double &t) const;


  virtual void postprocess_newly_created_triangulation(Triangulation<dim, spacedim> &tria) const;

  virtual void apply_dirichlet_bcs (const DoFHandler<dim,spacedim> &dof_handler,
                                    ConstraintMatrix &constraints) const;

  template<typename Number>
  void apply_neumann_bcs (const typename DoFHandler<dim,spacedim>::active_cell_iterator &cell,
                          Scratch &scratch,
                          CopySystem &data,
                          std::vector<Number> &local_residual) const;


  template<typename Number>
  void apply_forcing_terms (const typename DoFHandler<dim,spacedim>::active_cell_iterator &cell,
                            Scratch &scratch,
                            CopySystem &data,
                            std::vector<Number> &local_residual) const;

  virtual void initialize_data(const typename LAC::VectorType &solution,
                               const typename LAC::VectorType &solution_dot,
                               const double t,
                               const double alpha) const;


  virtual void get_preconditioner_energy(const typename DoFHandler<dim,spacedim>::active_cell_iterator &,
                                         Scratch &,
                                         CopySystem &,
                                         Sdouble &) const;

  virtual void get_preconditioner_energy(const typename DoFHandler<dim,spacedim>::active_cell_iterator &,
                                         Scratch &,
                                         CopyPreconditioner &,
                                         SSdouble &)  const;

  virtual void get_system_energy(const typename DoFHandler<dim,spacedim>::active_cell_iterator &,
                                 Scratch &,
                                 CopySystem &,
                                 Sdouble &) const;

  virtual void get_system_energy(const typename DoFHandler<dim,spacedim>::active_cell_iterator &,
                                 Scratch &,
                                 CopySystem &,
                                 SSdouble &)  const;

  virtual void get_system_residual (const typename DoFHandler<dim,spacedim>::active_cell_iterator &cell,
                                    Scratch &scratch,
                                    CopySystem &data,
                                    std::vector<Sdouble> &local_residual) const;

  virtual void get_system_residual (const typename DoFHandler<dim,spacedim>::active_cell_iterator &cell,
                                    Scratch &scratch,
                                    CopySystem &data,
                                    std::vector<double> &local_residual) const;

  virtual void get_preconditioner_residual (const typename DoFHandler<dim,spacedim>::active_cell_iterator &cell,
                                            Scratch &scratch,
                                            CopyPreconditioner &data,
                                            std::vector<Sdouble> &local_residual) const;

  virtual void compute_system_operators(const DoFHandler<dim,spacedim> &,
                                        const typename LAC::BlockMatrix &,
                                        const typename LAC::BlockMatrix &,
                                        const std::vector<shared_ptr<typename LAC::BlockMatrix> >,
                                        LinearOperator<typename LAC::VectorType> &,
                                        LinearOperator<typename LAC::VectorType> &) const;


  virtual void assemble_local_system (const typename DoFHandler<dim,spacedim>::active_cell_iterator &cell,
                                      Scratch &scratch,
                                      CopySystem &data) const;

  virtual void assemble_local_preconditioner (const typename DoFHandler<dim,spacedim>::active_cell_iterator &cell,
                                              Scratch &scratch,
                                              CopyPreconditioner &data) const;

  virtual shared_ptr<Mapping<dim,spacedim> > get_mapping(const DoFHandler<dim,spacedim> &,
                                                         const typename LAC::VectorType &) const;

  virtual UpdateFlags get_jacobian_flags() const;

  virtual UpdateFlags get_residual_flags() const;

  virtual UpdateFlags get_jacobian_preconditioner_flags() const;


  virtual UpdateFlags get_face_flags() const;

  void fix_solution_dot_derivative(FEValuesCache<dim,spacedim> &, double) const;

  void fix_solution_dot_derivative(FEValuesCache<dim,spacedim> &fe_cache, Sdouble alpha) const;


  void fix_solution_dot_derivative(FEValuesCache<dim,spacedim> &fe_cache, SSdouble alpha) const;

  template<typename Number>
  void reinit(const Number &alpha,
              const typename DoFHandler<dim,spacedim>::active_cell_iterator &cell,
              FEValuesCache<dim,spacedim> &fe_cache) const;


  template<typename Number>
  void reinit(const Number &alpha,
              const typename DoFHandler<dim,spacedim>::active_cell_iterator &cell,
              const unsigned int face_no,
              FEValuesCache<dim,spacedim> &fe_cache) const;


  // auxiliary matrices ///////////////////////////////////////////////////////
  // TODO aux_mtrices energy

  virtual void get_aux_matrix_residuals (const typename DoFHandler<dim,spacedim>::active_cell_iterator &cell,
                                         Scratch &scratch,
                                         CopyPreconditioner &data,
                                         std::vector<std::vector<double> > &local_residuals) const;

  virtual void get_aux_matrix_residuals (const typename DoFHandler<dim,spacedim>::active_cell_iterator &cell,
                                         Scratch &scratch,
                                         CopyPreconditioner &data,
                                         std::vector<std::vector<Sdouble> > &local_residuals) const;

  virtual void assemble_local_aux_matrices (const typename DoFHandler<dim,spacedim>::active_cell_iterator &cell,
                                            Scratch &scratch,
                                            CopyPreconditioner &data) const;

  //  virtual void set_aux_matrix_update_flags();

  UpdateFlags get_aux_matrix_flags(const unsigned int &i) const;

  const Table<2,DoFTools::Coupling> &get_aux_matrix_coupling(const unsigned int &i) const;
  virtual unsigned int get_number_of_aux_matrices() const;

  std::vector<UpdateFlags> aux_matrix_update_flags;
  std::vector<Table<2,DoFTools::Coupling> > aux_matrix_coupling;



protected:



  mutable ParsedMappedFunctions<spacedim,n_components>  forcing_terms; // on the volume
  mutable ParsedMappedFunctions<spacedim,n_components>  neumann_bcs;
  mutable ParsedDirichletBCs<dim,spacedim,n_components> dirichlet_bcs;

  std::string str_diff_comp;
  std::vector<unsigned int> _diff_comp;

  mutable const typename LAC::VectorType *solution;

  mutable typename LAC::VectorType old_solution;

  mutable const typename LAC::VectorType *solution_dot;

  mutable double t;

  mutable double old_t;

  mutable double alpha;
  mutable unsigned int dofs_per_cell;
  mutable unsigned int n_q_points;
  mutable unsigned int n_face_q_points;


};

template <int dim, int spacedim, int n_components, typename LAC>
template<typename Number>
void
Interface<dim,spacedim,n_components,LAC>::
apply_neumann_bcs (
  const typename DoFHandler<dim,spacedim>::active_cell_iterator &cell,
  Scratch &scratch,
  CopySystem &data,
  std::vector<Number> &local_residual) const
{

  Number dummy = 0.0;

  for (unsigned int face=0; face < GeometryInfo<dim>::faces_per_cell; ++face)
    {
      unsigned int face_id = cell->face(face)->boundary_id();
      if (cell->face(face)->at_boundary() && neumann_bcs.acts_on_id(face_id))
        {
          this->reinit(dummy, cell, face, scratch);

          auto &fev = scratch.get_current_fe_values();
          auto &q_points = scratch.get_quadrature_points();
          auto &JxW = scratch.get_JxW_values();

          for (unsigned int q=0; q<q_points.size(); ++q)
            {
              Vector<double> T(n_components);
              neumann_bcs.get_mapped_function(face_id)->vector_value(q_points[q], T);

              for (unsigned int i=0; i<local_residual.size(); ++i)
                for (unsigned int c=0; c<n_components; ++c)
                  local_residual[i] -= T[c]*fev.shape_value_component(i,q,c)*JxW[q];

            }
          break;
        }
    }
}

template <int dim, int spacedim, int n_components, typename LAC>
template<typename Number>
void
Interface<dim,spacedim,n_components,LAC>::
apply_forcing_terms (const typename DoFHandler<dim,spacedim>::active_cell_iterator &cell,
                     Scratch &scratch,
                     CopySystem &data,
                     std::vector<Number> &local_residual) const
{
  unsigned cell_id = cell->material_id();
  if (forcing_terms.acts_on_id(cell_id))
    {
      Number dummy = 0.0;
      this->reinit(dummy, cell, scratch);

      auto &fev = scratch.get_current_fe_values();
      auto &q_points = scratch.get_quadrature_points();
      auto &JxW = scratch.get_JxW_values();
      for (unsigned int q=0; q<q_points.size(); ++q)
        for (unsigned int i=0; i<local_residual.size(); ++i)
          for (unsigned int c=0; c<n_components; ++c)
            {
              double B = forcing_terms.get_mapped_function(cell_id)->value(q_points[q],c);
              local_residual[i] -= B*fev.shape_value_component(i,q,c)*JxW[q];
            }
    }
}


template <int dim, int spacedim, int n_components, typename LAC>
template<typename Number>
void
Interface<dim,spacedim,n_components,LAC>::reinit(const Number &alpha,
                                                 const typename DoFHandler<dim,spacedim>::active_cell_iterator &cell,
                                                 FEValuesCache<dim,spacedim> &fe_cache) const
{
  fe_cache.reinit(cell);
  fe_cache.cache_local_solution_vector("old_solution", this->old_solution, alpha);
  fe_cache.cache_local_solution_vector("solution", *this->solution, alpha);
  fe_cache.cache_local_solution_vector("solution_dot", *this->solution_dot, alpha);
  this->fix_solution_dot_derivative(fe_cache, alpha);
}


template <int dim, int spacedim, int n_components, typename LAC>
template<typename Number>
void
Interface<dim,spacedim,n_components,LAC>::reinit(const Number &alpha,
                                                 const typename DoFHandler<dim,spacedim>::active_cell_iterator &cell,
                                                 const unsigned int face_no,
                                                 FEValuesCache<dim,spacedim> &fe_cache) const
{
  fe_cache.reinit(cell, face_no);
  fe_cache.cache_local_solution_vector("old_solution", this->old_solution, alpha);
  fe_cache.cache_local_solution_vector("solution", *this->solution, alpha);
  fe_cache.cache_local_solution_vector("solution_dot", *this->solution_dot, alpha);
  this->fix_solution_dot_derivative(fe_cache, alpha);
}



#endif