Home
Publications
Presentations
Support
Download
Examples
Developers
Class Docs
Namespaces
Classes
Files
Mailing Lists
Gallery
Wiki

system.h

Go to the documentation of this file.

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

// Local Includes
#include "libmesh/auto_ptr.h"
#include "libmesh/elem_range.h"
#include "libmesh/enum_norm_type.h"
#include "libmesh/enum_xdr_mode.h"
#include "libmesh/enum_subset_solve_mode.h"
#include "libmesh/enum_parallel_type.h"
#include "libmesh/fe_type.h"
#include "libmesh/libmesh_common.h"
#include "libmesh/tensor_value.h" // For point_hessian
#include "libmesh/qoi_set.h"
#include "libmesh/reference_counted_object.h"
#include "libmesh/variable.h"
#include "libmesh/fem_function_base.h"

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

namespace libMesh
{

// Forward Declarations
class System;
class SystemNorm;
class EquationSystems;
class MeshBase;
class Xdr;
class DofMap;
template <typename Output> class FunctionBase;
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;
class SystemSubset;

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

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

  class Initialization
  {
  public:
    virtual ~Initialization () {}

    virtual void initialize () = 0;
  };



  class Assembly
  {
  public:
    virtual ~Assembly () {}

    virtual void assemble () = 0;
  };



  class Constraint
  {
  public:
    virtual ~Constraint () {}

    virtual void constrain () = 0;
  };



  class QOI
  {
  public:
    virtual ~QOI () {}

    virtual void qoi (const QoISet& qoi_indices) = 0;
  };



  class QOIDerivative
  {
  public:
    virtual ~QOIDerivative () {}

    virtual void qoi_derivative (const QoISet& qoi_indices) = 0;
  };



  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 restrict_solve_to (const SystemSubset* subset,
                                  const SubsetSolveMode subset_solve_mode=SUBSET_ZERO);

  virtual void solve () {}

  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());
  bool is_adjoint_already_solved() const
  { return adjoint_already_solved;}

  void set_adjoint_already_solved(bool setting)
  { adjoint_already_solved = setting;}


  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 "Basic"; }

  void project_solution (FunctionBase<Number> *f,
                         FunctionBase<Gradient> *g = NULL) const;

  void project_solution (FEMFunctionBase<Number> *f,
                         FEMFunctionBase<Gradient> *g = NULL) const;

  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),
                         const Parameters& parameters) const;

  void project_vector (NumericVector<Number>& new_vector,
                       FunctionBase<Number> *f,
                       FunctionBase<Gradient> *g = NULL) const;

  void project_vector (NumericVector<Number>& new_vector,
                       FEMFunctionBase<Number> *f,
                       FEMFunctionBase<Gradient> *g = NULL) 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),
                       const Parameters& parameters,
                       NumericVector<Number>& new_vector) const;

  void boundary_project_solution (const std::set<boundary_id_type> &b,
                                  const std::vector<unsigned int> &variables,
                                  FunctionBase<Number> *f,
                                  FunctionBase<Gradient> *g = NULL);

  void boundary_project_solution (const std::set<boundary_id_type> &b,
                                  const std::vector<unsigned int> &variables,
                                  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),
                                  const Parameters& parameters);

  void boundary_project_vector (const std::set<boundary_id_type> &b,
                                const std::vector<unsigned int> &variables,
                                NumericVector<Number>& new_vector,
                                FunctionBase<Number> *f,
                                FunctionBase<Gradient> *g = NULL) const;

  void boundary_project_vector (const std::set<boundary_id_type> &b,
                                const std::vector<unsigned int> &variables,
                                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),
                                const 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 ();

  void set_basic_system_only ();

  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,
                                      const ParallelType type = PARALLEL);

  void remove_vector(const std::string& vec_name);

  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) const;

  const std::string & vector_name (const NumericVector<Number> & vec_reference) const;

  void set_vector_preservation (const std::string &vec_name, bool preserve);

  bool vector_preservation (const std::string &vec_name) const;

  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;

  virtual unsigned int n_matrices () const;

  unsigned int n_vars() const;

  unsigned int n_variable_groups() const;

  unsigned int n_components() const;

  dof_id_type n_dofs() const;

  dof_id_type n_active_dofs() const;

  dof_id_type n_constrained_dofs() const;

  dof_id_type n_local_constrained_dofs() const;

  dof_id_type n_local_dofs() const;

  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);

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

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

  const Variable & variable (unsigned int var) const;

  const VariableGroup & variable_group (unsigned int vg) 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;

  void get_all_variable_numbers(std::vector<unsigned int>& all_variable_numbers) const;

  unsigned int variable_scalar_number (const std::string& var,
                                       unsigned int component) const;

  unsigned int variable_scalar_number (unsigned int var_num,
                                       unsigned int component) const;


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

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

  bool identify_variable_groups () const;

  void identify_variable_groups (const bool);

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

  Real calculate_norm(const 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);

  dof_id_type read_serialized_vectors (Xdr &io,
                                       const std::vector<NumericVector<Number>*> &vectors) const;
  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;

  dof_id_type write_serialized_vectors (Xdr &io,
                                        const std::vector<const NumericVector<Number>*> &vectors) 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_init_object (Initialization& init);

  void attach_assemble_function (void fptr(EquationSystems& es,
                                           const std::string& name));

  void attach_assemble_object (Assembly& assemble);

  void attach_constraint_function (void fptr(EquationSystems& es,
                                             const std::string& name));

  void attach_constraint_object (Constraint& constrain);

  void attach_QOI_function (void fptr(EquationSystems& es,
                                      const std::string& name,
                                      const QoISet& qoi_indices));

  void attach_QOI_object (QOI& qoi);

  void attach_QOI_derivative (void fptr(EquationSystems& es,
                                        const std::string& name,
                                        const QoISet& qoi_indices));

  void attach_QOI_derivative_object (QOIDerivative& qoi_derivative);

  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;

  virtual void disable_cache ();

  bool use_fixed_solution;

  int extra_quadrature_order;


  //--------------------------------------------------
  // The solution and solution access members

  Number current_solution (const dof_id_type global_dof_number) const;

  AutoPtr<NumericVector<Number> > solution;

  AutoPtr<NumericVector<Number> > current_local_solution;

  Real time;

  std::vector<Number> qoi;

  Number point_value(unsigned int var, const Point &p,
                     const bool insist_on_success = true) const;

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

  Gradient point_gradient(unsigned int var, const Point &p,
                          const bool insist_on_success = true) const;

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

  Tensor point_hessian(unsigned int var, const Point &p,
                       const bool insist_on_success = true) const;

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

  void local_dof_indices (const unsigned int var,
                          std::set<dof_id_type> & 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:
  System (const System&);

  System& operator=(const System&);

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

  template <typename iterator_type>
  dof_id_type read_serialized_blocked_dof_objects (const dof_id_type n_objects,
                                                   const iterator_type begin,
                                                   const iterator_type end,
                                                   Xdr &io,
                                                   const std::vector<NumericVector<Number>*> &vecs,
                                                   const unsigned int var_to_read=libMesh::invalid_uint) const;

  unsigned int read_SCALAR_dofs (const unsigned int var,
                                 Xdr &io,
                                 NumericVector<Number> &vec) const;

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

  template <typename iterator_type>
  dof_id_type write_serialized_blocked_dof_objects (const std::vector<const NumericVector<Number>*> &vecs,
                                                    const dof_id_type n_objects,
                                                    const iterator_type begin,
                                                    const iterator_type end,
                                                    Xdr &io,
                                                    const unsigned int var_to_write=libMesh::invalid_uint) const;

  unsigned int write_SCALAR_dofs (const NumericVector<Number> &vec,
                                  const unsigned int var,
                                  Xdr &io) const;

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

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

  Initialization * _init_system_object;

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

  Assembly * _assemble_system_object;

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

  Constraint * _constrain_system_object;

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

  QOI *_qoi_evaluate_object;

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

  QOIDerivative *_qoi_evaluate_derivative_object;

  AutoPtr<DofMap> _dof_map;

  EquationSystems& _equation_systems;

  MeshBase& _mesh;

  const std::string _sys_name;

  const unsigned int _sys_number;

  std::vector<Variable> _variables;

  std::vector<VariableGroup> _variable_groups;

  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;

  std::map<std::string, ParallelType> _vector_types;

  bool _solution_projection;

  bool _basic_system_only;

  bool _can_add_vectors;

  bool _identify_variable_groups;

  bool _additional_data_written;

  std::vector<unsigned int> _written_var_indices;

  bool adjoint_already_solved;
};



// ------------------------------------------------------------
// 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
void System::set_basic_system_only ()
{
  _basic_system_only = true;
}



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



inline
unsigned int System::n_variable_groups() const
{
  return libmesh_cast_int<unsigned int>(_variable_groups.size());
}



inline
unsigned int System::n_components() const
{
  if (_variables.empty())
    return 0;

  const Variable& last = _variables.back();
  return last.first_scalar_number() + last.n_components();
}



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

  return _variables[i];
}



inline
const VariableGroup & System::variable_group (const unsigned int vg) const
{
  libmesh_assert_less (vg, _variable_groups.size());

  return _variable_groups[vg];
}



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

  return _variables[i].name();
}



inline
unsigned int
System::variable_scalar_number (const std::string& var,
                                unsigned int component) const
{
  return variable_scalar_number(this->variable_number(var), component);
}



inline
unsigned int
System::variable_scalar_number (unsigned int var_num,
                                unsigned int component) const
{
  return _variables[var_num].first_scalar_number() + component;
}



inline
const FEType & System::variable_type (const unsigned int i) const
{
  libmesh_assert_less (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
bool System::identify_variable_groups () const
{
  return _identify_variable_groups;
}



inline
void System::identify_variable_groups (const bool ivg)
{
  _identify_variable_groups = ivg;
}



inline
dof_id_type 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 libmesh_cast_int<unsigned int>(_vectors.size());
}

inline
unsigned int System::n_matrices () const
{
  return 0;
}

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
void System::disable_cache () { assemble_before_solve = true; }

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();
}


} // namespace libMesh

#endif // LIBMESH_SYSTEM_H

---

Site Created By: libMesh Developers
Last modified: February 05 2013 19:54:48 UTC

Hosted By: