system.h

Go to the documentation of this file.

// $Id: system.h 3557 2009-11-12 00:05:53Z roystgnr $

// The libMesh Finite Element Library.
// Copyright (C) 2002-2008 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 __system_h__
#define __system_h__

// C++ includes
#include <vector>
#include <set>

// Local Includes
#include "auto_ptr.h"
#include "elem_range.h"
#include "enum_xdr_mode.h"
#include "fe_type.h"
#include "libmesh_common.h"
#include "qoi_set.h"
#include "reference_counted_object.h"
#include "system_norm.h"

// Forward Declarations
class System;
class EquationSystems;
class MeshBase;
class Xdr;
class DofMap;
class Parameters;
class ParameterVector;
class Point;
class SensitivityData;
template <typename T> class NumericVector;
template <typename T> class VectorValue;
typedef VectorValue<Number> NumberVectorValue;
typedef NumberVectorValue Gradient;

// ------------------------------------------------------------
// System class definition
class System : public ReferenceCountedObject<System>
{
protected:

  System (EquationSystems& es,
          const std::string& name,
          const unsigned int number);
  
public:

  
  virtual ~System ();

  typedef System sys_type;

  sys_type & system () { return *this; }
  
  virtual void clear ();
  
  void init ();
  
  virtual void reinit ();
  
  virtual void update ();
  
  virtual void assemble ();

  virtual void assemble_qoi
    (const QoISet &qoi_indices = QoISet());
 
  virtual void assemble_qoi_derivative 
    (const QoISet &qoi_indices = QoISet());
 
  virtual void assemble_residual_derivatives (const ParameterVector& parameters);
 
  virtual void solve () = 0;
  
  virtual std::pair<unsigned int, Real>
    sensitivity_solve (const ParameterVector& parameters);
  
  virtual std::pair<unsigned int, Real>
    weighted_sensitivity_solve (const ParameterVector& parameters,
                                const ParameterVector& weights);
 
  virtual std::pair<unsigned int, Real>
    adjoint_solve (const QoISet& qoi_indices = QoISet());
  
  virtual std::pair<unsigned int, Real>
    weighted_sensitivity_adjoint_solve (const ParameterVector& parameters,
                                        const ParameterVector& weights,
                                        const QoISet& qoi_indices = QoISet());
 
  virtual void qoi_parameter_sensitivity (const QoISet& qoi_indices,
                                          const ParameterVector& parameters,
                                          SensitivityData& sensitivities);
  
  virtual void adjoint_qoi_parameter_sensitivity (const QoISet& qoi_indices,
                                                  const ParameterVector& parameters,
                                                  SensitivityData& sensitivities);
  
  virtual void forward_qoi_parameter_sensitivity (const QoISet& qoi_indices,
                                                  const ParameterVector& parameters,
                                                  SensitivityData& sensitivities);
  
  virtual void qoi_parameter_hessian(const QoISet& qoi_indices,
                                     const ParameterVector& parameters,
                                     SensitivityData& hessian);

  virtual void qoi_parameter_hessian_vector_product(const QoISet& qoi_indices,
                                                    const ParameterVector& parameters,
                                                    const ParameterVector& vector,
                                                    SensitivityData& product);
  
  virtual bool compare (const System& other_system, 
                        const Real threshold,
                        const bool verbose) const;

  const std::string & name () const;

  virtual std::string system_type () const { return "BasicSystem"; }

  void project_solution (Number fptr(const Point& p,
                                     const Parameters& parameters,
                                     const std::string& sys_name,
                                     const std::string& unknown_name),
                         Gradient gptr(const Point& p,
                                       const Parameters& parameters,
                                       const std::string& sys_name,
                                       const std::string& unknown_name),
                         Parameters& parameters) const;

  void project_vector (Number fptr(const Point& p,
                                   const Parameters& parameters,
                                   const std::string& sys_name,
                                   const std::string& unknown_name),
                       Gradient gptr(const Point& p,
                                     const Parameters& parameters,
                                     const std::string& sys_name,
                                     const std::string& unknown_name),
                       Parameters& parameters,
                       NumericVector<Number>& new_vector) const;
  
  unsigned int number () const;
  
  void update_global_solution (std::vector<Number>& global_soln) const;

  void update_global_solution (std::vector<Number>& global_soln,
                               const unsigned int dest_proc) const;

  const MeshBase & get_mesh() const;
  
  MeshBase & get_mesh();
  
  const DofMap & get_dof_map() const;
  
  DofMap & get_dof_map();

  const EquationSystems & get_equation_systems() const { return _equation_systems; }

  EquationSystems & get_equation_systems() { return _equation_systems; }

  bool active () const;

  void activate ();

  void deactivate ();

  typedef std::map<std::string, NumericVector<Number>* >::iterator       vectors_iterator;
  typedef std::map<std::string, NumericVector<Number>* >::const_iterator const_vectors_iterator;

  vectors_iterator vectors_begin ();

  const_vectors_iterator vectors_begin () const;

  vectors_iterator vectors_end ();

  const_vectors_iterator vectors_end () const;

  NumericVector<Number> & add_vector (const std::string& vec_name,
                                      const bool projections=true);

  bool& project_solution_on_reinit (void)
    { return _solution_projection; }
  
  bool have_vector (const std::string& vec_name) const;
  
  const NumericVector<Number> * request_vector (const std::string& vec_name) const;
  
  NumericVector<Number> * request_vector (const std::string& vec_name);
  
  const NumericVector<Number> * request_vector (const unsigned int vec_num) const;

  NumericVector<Number> * request_vector (const unsigned int vec_num);

  const NumericVector<Number> & get_vector (const std::string& vec_name) const;

  NumericVector<Number> & get_vector (const std::string& vec_name);

  const NumericVector<Number> & get_vector (const unsigned int vec_num) const;

  NumericVector<Number> & get_vector (const unsigned int vec_num);

  const std::string & vector_name (const unsigned int vec_num);

  NumericVector<Number> & add_adjoint_solution(unsigned int i=0);

  NumericVector<Number> & get_adjoint_solution(unsigned int i=0);

  const NumericVector<Number> & get_adjoint_solution(unsigned int i=0) const;

  NumericVector<Number> & add_sensitivity_solution(unsigned int i=0);

  NumericVector<Number> & get_sensitivity_solution(unsigned int i=0);

  const NumericVector<Number> & get_sensitivity_solution(unsigned int i=0) const;

  NumericVector<Number> & add_weighted_sensitivity_adjoint_solution(unsigned int i=0);

  NumericVector<Number> & get_weighted_sensitivity_adjoint_solution(unsigned int i=0);

  const NumericVector<Number> & get_weighted_sensitivity_adjoint_solution(unsigned int i=0) const;

  NumericVector<Number> & add_weighted_sensitivity_solution();

  NumericVector<Number> & get_weighted_sensitivity_solution();

  const NumericVector<Number> & get_weighted_sensitivity_solution() const;

  NumericVector<Number> & add_adjoint_rhs(unsigned int i=0);

  NumericVector<Number> & get_adjoint_rhs(unsigned int i=0);

  const NumericVector<Number> & get_adjoint_rhs(unsigned int i=0) const;

  NumericVector<Number> & add_sensitivity_rhs(unsigned int i=0);

  NumericVector<Number> & get_sensitivity_rhs(unsigned int i=0);

  const NumericVector<Number> & get_sensitivity_rhs(unsigned int i=0) const;

  unsigned int n_vectors () const;

  unsigned int n_vars() const;

  unsigned int n_dofs() const;

  unsigned int n_active_dofs() const;

  unsigned int n_constrained_dofs() const;
  
  unsigned int n_local_dofs() const;

  class Variable
  {
  public:
    
    Variable (const std::string &var_name,
              const unsigned int var_number,
              const FEType &var_type) :
      _name(var_name),
      _number(var_number),
      _type(var_type)
    {}
    
    Variable (const std::string &var_name,
              const unsigned int var_number,
              const FEType &var_type,
              const std::set<subdomain_id_type> &var_active_subdomains) :
      _name(var_name),
      _number(var_number),
      _type(var_type),
      _active_subdomains(var_active_subdomains)
    {}
    
    const std::string & name() const 
    { return _name; }

    unsigned int number() const 
    { return _number; }

    const FEType & type() const 
    { return _type; }

    bool active_on_subdomain (const subdomain_id_type sid) const
    { return (_active_subdomains.empty() || _active_subdomains.count(sid));  }

    bool implicitly_active () const
    { return _active_subdomains.empty(); }
    
  private:
    std::string             _name; 
    unsigned int            _number;
    FEType                  _type;
    std::set<subdomain_id_type> _active_subdomains;
  };

  unsigned int add_variable (const std::string& var,
                             const FEType& type,
                             const std::set<subdomain_id_type> * const active_subdomains = NULL);

  unsigned int add_variable (const std::string& var,
                             const Order order = FIRST,
                             const FEFamily = LAGRANGE,
                             const std::set<subdomain_id_type> * const active_subdomains = NULL);

  const Variable & variable (unsigned int var) const;

  bool has_variable(const std::string& var) const;
  
  const std::string & variable_name(const unsigned int i) const;
  
  unsigned short int variable_number (const std::string& var) const;

  const FEType & variable_type (const unsigned int i) const;

  const FEType & variable_type (const std::string& var) const;

  Real calculate_norm(NumericVector<Number>& v,
                      unsigned int var = 0,
                      FEMNormType norm_type = L2) const;

  Real calculate_norm(NumericVector<Number>& v,
                      const SystemNorm &norm) const;

  void read_header (Xdr& io,
                    const std::string &version,
                    const bool read_header=true,
                    const bool read_additional_data=true,
                    const bool read_legacy_format=false);

  void read_legacy_data (Xdr& io,
                         const bool read_additional_data=true);

  void read_serialized_data (Xdr& io,
                             const bool read_additional_data=true);

  void read_parallel_data (Xdr &io,
                           const bool read_additional_data);
  void write_header (Xdr& io,
                     const std::string &version,
                     const bool write_additional_data) const;

  void write_serialized_data (Xdr& io,
                              const bool write_additional_data = true) const;

  void write_parallel_data (Xdr &io,
                            const bool write_additional_data) const;
  
  std::string get_info () const;

  void attach_init_function (void fptr(EquationSystems& es,
                                       const std::string& name));
  
  void attach_assemble_function (void fptr(EquationSystems& es,
                                           const std::string& name));
  
  void attach_constraint_function (void fptr(EquationSystems& es,
                                             const std::string& name));
  
  void attach_QOI_function (void fptr(EquationSystems& es,
                                      const std::string& name,
                                      const QoISet& qoi_indices));
  
  void attach_QOI_derivative (void fptr(EquationSystems& es,
                                        const std::string& name,
                                        const QoISet& qoi_indices));
  
  virtual void user_initialization ();
  
  virtual void user_assembly ();
  
  virtual void user_constrain ();

  virtual void user_QOI (const QoISet& qoi_indices);

  virtual void user_QOI_derivative (const QoISet& qoi_indices);

  virtual void re_update ();

  virtual void restrict_vectors ();

  virtual void prolong_vectors ();

  bool assemble_before_solve;


  //--------------------------------------------------
  // The solution and solution access members
  
  Number current_solution (const unsigned int global_dof_number) const;

  AutoPtr<NumericVector<Number> > solution;

  AutoPtr<NumericVector<Number> > current_local_solution;

  std::vector<Number> qoi;

  void local_dof_indices (const unsigned int var,
                          std::set<unsigned int> & var_indices) const;

  void zero_variable (NumericVector<Number>& v, unsigned int var_num) const;

protected:

  virtual void init_data ();

  void project_vector (NumericVector<Number>&) const;

  void project_vector (const NumericVector<Number>&,
                       NumericVector<Number>&) const;
  
private:

  Real discrete_var_norm(NumericVector<Number>& v,
                         unsigned int var,
                         FEMNormType norm_type) const;

  template <typename iterator_type>
  unsigned int read_serialized_blocked_dof_objects (const unsigned int var,
                                                    const unsigned int n_objects,
                                                    const iterator_type begin,
                                                    const iterator_type end,
                                                    Xdr &io,
                                                    NumericVector<Number> &vec) const;

  void read_serialized_vector (Xdr& io,
                               NumericVector<Number> &vec);

  template <typename iterator_type>
  unsigned int write_serialized_blocked_dof_objects (const NumericVector<Number> &vec,
                                                     const unsigned int var,
                                                     const unsigned int n_objects,
                                                     const iterator_type begin,
                                                     const iterator_type end,
                                                     Xdr &io) const;

  void write_serialized_vector (Xdr& io,
                                const NumericVector<Number> &vec) const;

  void (* _init_system) (EquationSystems& es,
                         const std::string& name);
  
  void (* _assemble_system) (EquationSystems& es,
                             const std::string& name);

  void (* _constrain_system) (EquationSystems& es, 
                              const std::string& name);

  void (* _qoi_evaluate) (EquationSystems& es, 
                          const std::string& name,
                          const QoISet& qoi_indices);

  void (* _qoi_evaluate_derivative) (EquationSystems& es, 
                                     const std::string& name,
                                     const QoISet& qoi_indices);

  AutoPtr<DofMap> _dof_map;

  EquationSystems& _equation_systems;
  
  MeshBase& _mesh;
  
  const std::string _sys_name;

  const unsigned int _sys_number;
  
  std::vector<System::Variable> _variables;

  std::map<std::string, unsigned short int> _variable_numbers;

  bool _active;
  
  std::map<std::string, NumericVector<Number>* > _vectors;

  std::map<std::string, bool> _vector_projections;

  bool _solution_projection;

  bool _can_add_vectors;

  bool _additional_data_written;


  class ProjectVector 
  {
  private:
    const System                &system;
    const NumericVector<Number> &old_vector;
    NumericVector<Number>       &new_vector;

  public:
    ProjectVector (const System &system_in,
                   const NumericVector<Number> &old_v_in,
                   NumericVector<Number> &new_v_in) :
    system(system_in),
    old_vector(old_v_in),
    new_vector(new_v_in)
    {}
    
    void operator()(const ConstElemRange &range) const;
  };


  class BuildProjectionList
  {
  private:
    const System              &system;

  public:
    BuildProjectionList (const System &system_in) :
    system(system_in)
    {}

    BuildProjectionList (BuildProjectionList &other, Threads::split) :
     system(other.system)
    {}
    
    void unique();
    void operator()(const ConstElemRange &range);
    void join (const BuildProjectionList &other);
    std::vector<unsigned int> send_list;
  };



  class ProjectSolution
  {
  private:
    const System                &system;

    Number (* fptr)(const Point& p,
                    const Parameters& parameters,
                    const std::string& sys_name,
                    const std::string& unknown_name);

    Gradient (* gptr)(const Point& p,
                      const Parameters& parameters,
                      const std::string& sys_name,
                      const std::string& unknown_name);

    Parameters &parameters;
    NumericVector<Number>       &new_vector;

  public:
    ProjectSolution (const System &system_in,
                     Number fptr_in(const Point& p,
                                    const Parameters& parameters,
                                    const std::string& sys_name,
                                    const std::string& unknown_name),
                     Gradient gptr_in(const Point& p,
                                      const Parameters& parameters,
                                      const std::string& sys_name,
                                      const std::string& unknown_name),
                     Parameters &parameters_in,
                     NumericVector<Number> &new_v_in) :
    system(system_in),
    fptr(fptr_in),
    gptr(gptr_in),
    parameters(parameters_in),
    new_vector(new_v_in)
    {}
    
    void operator()(const ConstElemRange &range) const;
  };
};



// ------------------------------------------------------------
// System inline methods
inline
const std::string & System::name() const
{
  return _sys_name;
}



inline
unsigned int System::number() const
{
  return _sys_number;
}



inline
const MeshBase & System::get_mesh() const
{
  return _mesh;
}



inline
MeshBase & System::get_mesh()
{
  return _mesh;
}



inline
const DofMap & System::get_dof_map() const
{
  return *_dof_map;
}



inline
DofMap & System::get_dof_map()
{
  return *_dof_map;
}



inline
bool System::active() const
{
  return _active;  
}



inline
void System::activate ()
{
  _active = true;
}



inline
void System::deactivate ()
{
  _active = false;
}



inline
unsigned int System::n_vars() const
{
  return _variables.size();
}



inline
const System::Variable & System::variable (const unsigned int i) const
{
  libmesh_assert (i < _variables.size());

  return _variables[i];
}



inline
const std::string & System::variable_name (const unsigned int i) const
{
  libmesh_assert (i < _variables.size());

  return _variables[i].name();
}



inline
const FEType & System::variable_type (const unsigned int i) const
{
  libmesh_assert (i < _variables.size());
  
  return _variables[i].type();
}



inline
const FEType & System::variable_type (const std::string& var) const
{
  return _variables[this->variable_number(var)].type();
}



inline
unsigned int System::n_active_dofs() const
{
  return this->n_dofs() - this->n_constrained_dofs();
}



inline
bool System::have_vector (const std::string& vec_name) const
{
  return (_vectors.count(vec_name));
}



inline
unsigned int System::n_vectors () const
{
  return _vectors.size(); 
}

inline
System::vectors_iterator System::vectors_begin ()
{
  return _vectors.begin();
}

inline
System::const_vectors_iterator System::vectors_begin () const
{
  return _vectors.begin();
}

inline
System::vectors_iterator System::vectors_end ()
{
  return _vectors.end();
}

inline
System::const_vectors_iterator System::vectors_end () const
{
  return _vectors.end();
}

inline
void System::assemble_residual_derivatives (const ParameterVector&)
{
  libmesh_not_implemented();
}

inline
std::pair<unsigned int, Real>
System::sensitivity_solve (const ParameterVector&)
{
  libmesh_not_implemented();
}

inline
std::pair<unsigned int, Real>
System::weighted_sensitivity_solve (const ParameterVector&,
                                    const ParameterVector&)
{
  libmesh_not_implemented();
}

inline
std::pair<unsigned int, Real>
System::adjoint_solve (const QoISet&)
{
  libmesh_not_implemented();
}

inline
std::pair<unsigned int, Real>
System::weighted_sensitivity_adjoint_solve (const ParameterVector&,
                                            const ParameterVector&,
                                            const QoISet&)
{
  libmesh_not_implemented();
}

inline
void
System::adjoint_qoi_parameter_sensitivity (const QoISet&,
                                           const ParameterVector&,
                                           SensitivityData&)
{
  libmesh_not_implemented();
}

inline
void
System::forward_qoi_parameter_sensitivity (const QoISet&,
                                           const ParameterVector&,
                                           SensitivityData&)
{
  libmesh_not_implemented();
}

inline
void
System::qoi_parameter_hessian(const QoISet&,
                              const ParameterVector&,
                              SensitivityData&)
{
  libmesh_not_implemented();
}

inline
void
System::qoi_parameter_hessian_vector_product(const QoISet&,
                                             const ParameterVector&,
                                             const ParameterVector&,
                                             SensitivityData&)
{
  libmesh_not_implemented();
}

#endif // #define __system_h__