Poster of Linux kernelThe best gift for a Linux geek
ContinuationSystem

ContinuationSystem

Section: C Library Functions (3) Updated: Thu Apr 7 2011
Local index Up
 

NAME

ContinuationSystem -  

SYNOPSIS


#include <continuation_system.h>

Inherits FEMSystem.  

Public Types


enum Predictor { Euler, AB2, Invalid_Predictor }

typedef ContinuationSystem sys_type

typedef FEMSystem Parent

typedef std::map< std::string, SparseMatrix< Number > * >::iterator matrices_iterator

typedef std::map< std::string, SparseMatrix< Number > * >::const_iterator const_matrices_iterator

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

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

Public Member Functions


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

virtual ~ContinuationSystem ()

virtual void clear ()

virtual void solve ()

void continuation_solve ()

void advance_arcstep ()

void set_max_arclength_stepsize (Real maxds)

void save_current_solution ()

virtual void assembly (bool get_residual, bool get_jacobian)

virtual void time_evolving (unsigned int var)

virtual void mesh_x_position (unsigned int sysnum, unsigned int var)

virtual void mesh_y_position (unsigned int sysnum, unsigned int var)

virtual void mesh_z_position (unsigned int sysnum, unsigned int var)

void mesh_position_get ()

void mesh_position_set ()

virtual bool eulerian_residual (bool request_jacobian, DiffContext &context)

virtual bool mass_residual (bool request_jacobian, DiffContext &context)

virtual AutoPtr< DiffContext > build_context ()

virtual void init_context (DiffContext &)

virtual void postprocess ()

virtual void assemble_qoi ()

virtual void assemble_qoi_derivative ()

virtual void reinit ()

virtual void assemble ()

virtual bool element_time_derivative (bool request_jacobian, DiffContext &)

virtual bool element_constraint (bool request_jacobian, DiffContext &)

virtual bool side_time_derivative (bool request_jacobian, DiffContext &)

virtual bool side_constraint (bool request_jacobian, DiffContext &)

virtual void element_postprocess (DiffContext &)

virtual void side_postprocess (DiffContext &)

virtual void element_qoi (DiffContext &)

virtual void element_qoi_derivative (DiffContext &)

virtual void side_qoi (DiffContext &)

virtual void side_qoi_derivative (DiffContext &)

virtual bool side_mass_residual (bool request_jacobian, DiffContext &)

virtual void adjoint_solve ()

virtual void qoi_parameter_sensitivity (std::vector< Number * > &parameters, std::vector< Number > &sensitivities)

sys_type & system ()

virtual std::string system_type () const

SparseMatrix< Number > & add_matrix (const std::string &mat_name)

bool have_matrix (const std::string &mat_name) const

const SparseMatrix< Number > & get_matrix (const std::string &mat_name) const

SparseMatrix< Number > & get_matrix (const std::string &mat_name)

unsigned int n_matrices () const

void init ()

virtual void update ()

virtual bool compare (const System &other_system, const Real threshold, const bool verbose) const

const std::string & name () 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), 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

EquationSystems & get_equation_systems ()

bool active () const

void activate ()

void deactivate ()

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)

bool have_vector (const std::string &vec_name) const

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 > & get_adjoint_solution ()

const NumericVector< Number > & get_adjoint_solution () 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

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

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

virtual void user_initialization ()

virtual void user_assembly ()

virtual void user_constrain ()

virtual void user_QOI ()

virtual void user_QOI_derivative ()

virtual void re_update ()

virtual void restrict_vectors ()

virtual void prolong_vectors ()

Number current_solution (const unsigned int global_dof_number) const

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
 

Static Public Member Functions


static std::string get_info ()

static void print_info ()

static unsigned int n_objects ()
 

Public Attributes


Real * continuation_parameter

bool quiet

Real continuation_parameter_tolerance

Real solution_tolerance

Real initial_newton_tolerance

Real old_continuation_parameter

Real min_continuation_parameter

Real max_continuation_parameter

Real Theta

Real Theta_LOCA

unsigned int n_backtrack_steps

unsigned int n_arclength_reductions

Real ds_min

Predictor predictor

Real newton_stepgrowth_aggressiveness

bool newton_progress_check

bool fe_reinit_during_postprocess

int extra_quadrature_order

Real numerical_jacobian_h

Real verify_analytic_jacobians

bool compute_internal_sides

bool postprocess_sides

bool assemble_qoi_sides

AutoPtr< TimeSolver > time_solver

Real time

Real deltat

bool print_solution_norms

bool print_solutions

bool print_residual_norms

bool print_residuals

bool print_jacobian_norms

bool print_jacobians

bool print_element_jacobians

bool use_fixed_solution

SparseMatrix< Number > * matrix

NumericVector< Number > * rhs

bool assemble_before_solve

AutoPtr< NumericVector< Number > > solution

AutoPtr< NumericVector< Number > > current_local_solution

Number qoi
 

Protected Types


enum RHS_Mode { Residual, G_Lambda }

typedef bool(TimeSolver::* TimeSolverResPtr )(bool, DiffContext &)

typedef std::map< std::string, std::pair< unsigned int, unsigned int > > Counts
 

Protected Member Functions


virtual void init_data ()

void numerical_jacobian (TimeSolverResPtr res, FEMContext &context)

void numerical_elem_jacobian (FEMContext &context)

void numerical_side_jacobian (FEMContext &context)

virtual void init_matrices ()

void project_vector (NumericVector< Number > &) const

void project_vector (const NumericVector< Number > &, NumericVector< Number > &) const

void increment_constructor_count (const std::string &name)

void increment_destructor_count (const std::string &name)
 

Protected Attributes


RHS_Mode rhs_mode

unsigned int _mesh_sys

unsigned int _mesh_x_var

unsigned int _mesh_y_var

unsigned int _mesh_z_var

std::vector< bool > _time_evolving
 

Static Protected Attributes


static Counts _counts

static Threads::atomic< unsigned int > _n_objects

static Threads::spin_mutex _mutex
 

Private Member Functions


void initialize_tangent ()

void solve_tangent ()

void update_solution ()

void set_Theta ()

void set_Theta_LOCA ()

void apply_predictor ()
 

Private Attributes


NumericVector< Number > * du_ds

NumericVector< Number > * previous_du_ds

NumericVector< Number > * previous_u

NumericVector< Number > * y

NumericVector< Number > * y_old

NumericVector< Number > * z

NumericVector< Number > * delta_u

AutoPtr< LinearSolver< Number > > linear_solver

bool tangent_initialized

NewtonSolver * newton_solver

Real dlambda_ds

Real ds

Real ds_current

Real previous_dlambda_ds

Real previous_ds

unsigned int newton_step
 

Detailed Description

This class inherits from the FEMSystem. It can be used to do arclength continuation. Most of the ideas and the notation here come from HB Keller's 1977 paper:

{Kell-1977, author = {H.~B.~Keller}, title = {{Numerical solution of bifurcation and nonlinear eigenvalue problems}}, booktitle = {Applications of Bifurcation Theory, P.~H.~Rabinowitz (ed.)}, year = 1977, publisher = {Academic Press}, pages = {359--389}, notes = {QA 3 U45 No. 38 (PMA)} }

Author:

John W. Peterson 2007

Definition at line 55 of file continuation_system.h.  

Member Typedef Documentation

 

typedef std::map<std::string, SparseMatrix<Number>* >::const_iterator ImplicitSystem::const_matrices_iterator [inherited]

Definition at line 112 of file implicit_system.h.  

typedef std::map<std::string, NumericVector<Number>* >::const_iterator System::const_vectors_iterator [inherited]

Definition at line 279 of file system.h.  

typedef std::map<std::string, std::pair<unsigned int, unsigned int> > ReferenceCounter::Counts [protected, inherited]Data structure to log the information. The log is identified by the class name.

Definition at line 105 of file reference_counter.h.  

typedef std::map<std::string, SparseMatrix<Number>* >::iterator ImplicitSystem::matrices_iterator [inherited]Matrix iterator typedefs.

Definition at line 111 of file implicit_system.h.  

typedef FEMSystem ContinuationSystem::ParentThe type of the parent.

Reimplemented from FEMSystem.

Definition at line 79 of file continuation_system.h.  

typedef ContinuationSystem ContinuationSystem::sys_typeThe type of system.

Reimplemented from FEMSystem.

Definition at line 74 of file continuation_system.h.  

typedef bool(TimeSolver::* FEMSystem::TimeSolverResPtr)(bool, DiffContext &) [protected, inherited]Syntax sugar to make numerical_jacobian() declaration easier.

Definition at line 281 of file fem_system.h.  

typedef std::map<std::string, NumericVector<Number>* >::iterator System::vectors_iterator [inherited]Vector iterator typedefs.

Definition at line 278 of file system.h.  

Member Enumeration Documentation

 

enum ContinuationSystem::PredictorThe code provides the ability to select from different predictor schemes for getting the initial guess for the solution at the next point on the solution arc.

Enumerator:

Euler
First-order Euler predictor
AB2
Second-order explicit Adams-Bashforth predictor
Invalid_Predictor
Invalid predictor

Definition at line 221 of file continuation_system.h.

                 {
    Euler,

    AB2,

    Invalid_Predictor
  };
 

enum ContinuationSystem::RHS_Mode [protected]There are (so far) two different vectors which may be assembled using the assembly routine: 1.) The Residual = the normal PDE weighted residual 2.) G_Lambda = the derivative wrt the parameter lambda of the PDE weighted residual

It is up to the derived class to handle writing separate assembly code for the different cases. Usually something like: switch (rhs_mode) { case Residual: { Fu(i) += ... // normal PDE residual break; }

case G_Lambda: { Fu(i) += ... // derivative wrt control parameter break; }

Enumerator:

Residual
G_Lambda

Definition at line 286 of file continuation_system.h.

                {Residual,
                 G_Lambda};
 

Constructor & Destructor Documentation

 

ContinuationSystem::ContinuationSystem (EquationSystems &es, const std::string &name, const unsigned intnumber)Constructor. Optionally initializes required data structures.

Definition at line 27 of file continuation_system.C.

  : Parent(es, name, number),
    continuation_parameter(NULL),
    quiet(true),
    continuation_parameter_tolerance(1.e-6),
    solution_tolerance(1.e-6),
    initial_newton_tolerance(0.01),
    old_continuation_parameter(0.),
    min_continuation_parameter(0.),
    max_continuation_parameter(0.),
    Theta(1.),
    Theta_LOCA(1.),
    //tau(1.),
    n_backtrack_steps(5),
    n_arclength_reductions(5),
    ds_min(1.e-8),
    predictor(Euler),
    newton_stepgrowth_aggressiveness(1.),
    newton_progress_check(true),
    rhs_mode(Residual),
    linear_solver(LinearSolver<Number>::build()),
    tangent_initialized(false),
    newton_solver(NULL),
    dlambda_ds(0.707),
    ds(0.1),
    ds_current(0.1),
    previous_dlambda_ds(0.),
    previous_ds(0.),
    newton_step(0)
{
  // Warn about using untested code
  libmesh_experimental();
}
 

ContinuationSystem::~ContinuationSystem () [virtual]Destructor.

Definition at line 66 of file continuation_system.C.

References clear().

{
  this->clear();
}
 

Member Function Documentation

 

void System::activate () [inline, inherited]Activates the system. Only active systems are solved.

Definition at line 1103 of file system.h.

References System::_active.

{
  _active = true;
}
 

bool System::active () const [inline, inherited]Returns:

true if the system is active, false otherwise. An active system will be solved.

Definition at line 1095 of file system.h.

References System::_active.

{
  return _active;  
}
 

SparseMatrix< Number > & ImplicitSystem::add_matrix (const std::string &mat_name) [inherited]Adds the additional matrix mat_name to this system. Only allowed prior to assemble(). All additional matrices have the same sparsity pattern as the matrix used during solution. When not System but the user wants to initialize the mayor matrix, then all the additional matrices, if existent, have to be initialized by the user, too.

Definition at line 168 of file implicit_system.C.

References ImplicitSystem::_can_add_matrices, ImplicitSystem::_matrices, and ImplicitSystem::have_matrix().

Referenced by ImplicitSystem::add_system_matrix(), EigenTimeSolver::init(), and NewmarkSystem::NewmarkSystem().

{
  // only add matrices before initializing...
  if (!_can_add_matrices)
    {
      std::cerr << 'ERROR: Too late.  Cannot add matrices to the system after initialization'
                << std::endl
                << ' any more.  You should have done this earlier.'
                << std::endl;
      libmesh_error();
    }

  // Return the matrix if it is already there.
  if (this->have_matrix(mat_name))
    return *(_matrices[mat_name]);

  // Otherwise build the matrix and return it.
  SparseMatrix<Number>* buf = SparseMatrix<Number>::build().release();
  _matrices.insert (std::make_pair (mat_name, buf));

  return *buf;
}
 

unsigned int System::add_variable (const std::string &var, const FEType &type, const std::set< subdomain_id_type > *constactive_subdomains = NULL) [inherited]Adds the variable var to the list of variables for this system. Returns the index number for the new variable.

Definition at line 652 of file system.C.

References System::_dof_map, System::_variable_numbers, System::_variables, System::n_vars(), System::number(), System::variable_name(), and System::variable_type().

Referenced by System::add_variable(), ErrorVector::plot_error(), and System::read_header().

{  
  // Make sure the variable isn't there already
  // or if it is, that it's the type we want
  for (unsigned int v=0; v<this->n_vars(); v++)
    if (this->variable_name(v) == var)
      {
        if (this->variable_type(v) == type)
          return _variables[v].number();

        std::cerr << 'ERROR: incompatible variable '
                  << var
                  << ' has already been added for this system!'
                  << std::endl;
        libmesh_error();
      }

  const unsigned int curr_n_vars = this->n_vars();                                              

  // Add the variable to the list
  _variables.push_back((active_subdomains == NULL) ?
                       Variable(var, curr_n_vars, type) :
                       Variable(var, curr_n_vars, type, *active_subdomains));

  libmesh_assert ((curr_n_vars+1) == this->n_vars());

  _variable_numbers[var] = curr_n_vars;

  // Add the variable to the _dof_map
  _dof_map->add_variable (_variables.back());

  // Return the number of the new variable
  return curr_n_vars;
}
 

unsigned int System::add_variable (const std::string &var, const Orderorder = FIRST, const FEFamilyfamily = LAGRANGE, const std::set< subdomain_id_type > *constactive_subdomains = NULL) [inherited]Adds the variable var to the list of variables for this system. Same as before, but assumes LAGRANGE as default value for FEType.family.

Definition at line 691 of file system.C.

References System::add_variable().

{
  return this->add_variable(var, 
                            FEType(order, family), 
                            active_subdomains);
}
 

NumericVector< Number > & System::add_vector (const std::string &vec_name, const boolprojections = true) [inherited]Adds the additional vector vec_name to this system. Only allowed prior to init(). All the additional vectors are similarly distributed, like the solution, and inititialized to zero.

By default vectors added by add_vector are projected to changed grids by reinit(). To zero them instead (more efficient), pass 'false' as the second argument

Definition at line 530 of file system.C.

References System::_can_add_vectors, System::_vector_projections, System::_vectors, System::have_vector(), NumericVector< T >::init(), System::n_dofs(), System::n_local_dofs(), and libMeshEnums::PARALLEL.

Referenced by ExplicitSystem::add_system_rhs(), NonlinearImplicitSystem::adjoint_solve(), NewtonSolver::adjoint_solve(), LinearImplicitSystem::adjoint_solve(), UnsteadySolver::init(), init_data(), NewmarkSystem::NewmarkSystem(), System::read_header(), FrequencySystem::set_frequencies(), FrequencySystem::set_frequencies_by_range(), and FrequencySystem::set_frequencies_by_steps().

{
  // Return the vector if it is already there.
  if (this->have_vector(vec_name))
    return *(_vectors[vec_name]);

  // Otherwise build the vector
  NumericVector<Number>* buf = NumericVector<Number>::build().release();
  _vectors.insert (std::make_pair (vec_name, buf));
  _vector_projections.insert (std::make_pair (vec_name, projections));

  // Initialize it if necessary
  if (!_can_add_vectors)
    buf->init (this->n_dofs(), this->n_local_dofs(), false, PARALLEL);

  return *buf;
}
 

void DifferentiableSystem::adjoint_solve () [virtual, inherited]Invokes the adjoint solver associated with the system. For steady state solvers, this will find x where dF/dx = q. For transient solvers, this should integrate -dx/dt = dF/dx - q.

FIXME - transient adjoint solves are not yet implemented.

Reimplemented from ExplicitSystem.

Definition at line 112 of file diff_system.C.

References DifferentiableSystem::time_solver.

{
  time_solver->adjoint_solve();
}
 

void ContinuationSystem::advance_arcstep ()Call this function after a continuation solve to compute the tangent and get the next initial guess.

Definition at line 945 of file continuation_system.C.

References solve_tangent(), and update_solution().

{
  // Solve for the updated tangent du1/ds, d(lambda1)/ds
  solve_tangent();

  // Advance the solution and the parameter to the next value.
  update_solution();
}
 

void ContinuationSystem::apply_predictor () [private]Applies the predictor to the current solution to get a guess for the next solution.

Definition at line 1384 of file continuation_system.C.

References AB2, continuation_parameter, dlambda_ds, ds_current, du_ds, Euler, predictor, previous_dlambda_ds, previous_ds, previous_du_ds, and System::solution.

Referenced by continuation_solve(), and update_solution().

{
  if (predictor == Euler)
    {
      // 1.) Euler Predictor
      // Predict next the solution
      solution->add(ds_current, *du_ds);
      solution->close();
  
      // Predict next parameter value
      *continuation_parameter += ds_current*dlambda_ds;
    }


  else if (predictor == AB2)
    {
      // 2.) 2nd-order explicit AB predictor
      libmesh_assert(previous_ds != 0.0);
      const Real stepratio = ds_current/previous_ds;
      
      // Build up next solution value.
      solution->add( 0.5*ds_current*(2.+stepratio), *du_ds);
      solution->add(-0.5*ds_current*stepratio     , *previous_du_ds);
      solution->close();
      
      // Next parameter value
      *continuation_parameter +=
        0.5*ds_current*((2.+stepratio)*dlambda_ds -
                        stepratio*previous_dlambda_ds);
    }

  else
    {
      // Unknown predictor
      libmesh_error();
    }
  
}
 

void DifferentiableSystem::assemble () [virtual, inherited]Prepares matrix and rhs for matrix assembly. Users should not reimplement this

Reimplemented from ImplicitSystem.

Definition at line 98 of file diff_system.C.

References DifferentiableSystem::assembly().

{
  this->assembly(true, true);
}
 

void FEMSystem::assemble_qoi () [virtual, inherited]Runs a qoi assembly loop over all elements, and if assemble_qoi_sides is true over all sides.

Users may have to override this function for quantities of interest that are not expressible as a sum of element qois.

Reimplemented from ExplicitSystem.

Definition at line 486 of file fem_system.C.

References MeshBase::active_local_elements_begin(), MeshBase::active_local_elements_end(), DifferentiableSystem::assemble_qoi_sides, FEMSystem::build_context(), DifferentiableSystem::compute_internal_sides, FEMContext::elem, DifferentiableSystem::element_qoi(), System::get_mesh(), Elem::n_sides(), Elem::neighbor(), System::qoi, FEMContext::reinit(), FEMContext::side, FEMContext::side_fe, DifferentiableSystem::side_qoi(), and System::update().

{
  START_LOG('assemble_qoi()', 'FEMSystem');

  const MeshBase& mesh = this->get_mesh();

  this->update();

  AutoPtr<DiffContext> con = this->build_context();
  FEMContext &_femcontext = libmesh_cast_ref<FEMContext&>(*con);

  // the quantity of interest is assumed to be a sum of element and
  // side terms
  qoi = 0;

  // Loop over every active mesh element on this processor
  MeshBase::const_element_iterator el =
    mesh.active_local_elements_begin();
  const MeshBase::const_element_iterator end_el =
    mesh.active_local_elements_end();

  for ( ; el != end_el; ++el)
    {
      _femcontext.reinit(*this, *el);

      this->element_qoi(_femcontext);

      for (_femcontext.side = 0;
           _femcontext.side != _femcontext.elem->n_sides();
           ++_femcontext.side)
        {
          // Don't compute on non-boundary sides unless requested
          if (!assemble_qoi_sides ||
              (!compute_internal_sides &&
               _femcontext.elem->neighbor(_femcontext.side) != NULL))
            continue;

          std::map<FEType, FEBase *>::iterator fe_end =
            _femcontext.side_fe.end();
          for (std::map<FEType, FEBase *>::iterator i =
            _femcontext.side_fe.begin();
               i != fe_end; ++i)
            {
              i->second->reinit(_femcontext.elem, _femcontext.side);
            }

          this->side_qoi(_femcontext);
        }
    }

  Parallel::sum(qoi);

  STOP_LOG('assemble_qoi()', 'FEMSystem');
}
 

void FEMSystem::assemble_qoi_derivative () [virtual, inherited]Runs a qoi derivative assembly loop over all elements, and if assemble_qoi_sides is true over all sides.

Users may have to override this function for quantities of interest that are not expressible as a sum of element qois.

Reimplemented from ExplicitSystem.

Definition at line 543 of file fem_system.C.

References MeshBase::active_local_elements_begin(), MeshBase::active_local_elements_end(), NumericVector< T >::add_vector(), DifferentiableSystem::assemble_qoi_sides, FEMSystem::build_context(), DifferentiableSystem::compute_internal_sides, DofMap::constrain_element_vector(), DiffContext::dof_indices, FEMContext::elem, DiffContext::elem_residual, DifferentiableSystem::element_qoi_derivative(), System::get_dof_map(), System::get_mesh(), Elem::n_sides(), Elem::neighbor(), FEMContext::reinit(), ExplicitSystem::rhs, FEMContext::side, FEMContext::side_fe, DifferentiableSystem::side_qoi_derivative(), System::update(), and NumericVector< T >::zero().

{
  START_LOG('assemble_qoi_derivative()', 'FEMSystem');

  const MeshBase& mesh = this->get_mesh();

  this->update();

  AutoPtr<DiffContext> con = this->build_context();
  FEMContext &_femcontext = libmesh_cast_ref<FEMContext&>(*con);

  // In case there's already a rhs in use
  rhs->zero();

  // Loop over every active mesh element on this processor
  MeshBase::const_element_iterator el =
    mesh.active_local_elements_begin();
  const MeshBase::const_element_iterator end_el =
    mesh.active_local_elements_end();

  for ( ; el != end_el; ++el)
    {
      _femcontext.reinit(*this, *el);

      this->element_qoi_derivative(_femcontext);

      for (_femcontext.side = 0;
           _femcontext.side != _femcontext.elem->n_sides();
           ++_femcontext.side)
        {
          // Don't compute on non-boundary sides unless requested
          if (!assemble_qoi_sides ||
              (!compute_internal_sides &&
               _femcontext.elem->neighbor(_femcontext.side) != NULL))
            continue;

          std::map<FEType, FEBase *>::iterator fe_end =
            _femcontext.side_fe.end();
          for (std::map<FEType, FEBase *>::iterator i =
            _femcontext.side_fe.begin();
               i != fe_end; ++i)
            {
              i->second->reinit(_femcontext.elem, _femcontext.side);
            }

          this->side_qoi_derivative(_femcontext);
        }

      this->get_dof_map().constrain_element_vector
        (_femcontext.elem_residual, _femcontext.dof_indices, false);

      this->rhs->add_vector (_femcontext.elem_residual,
                             _femcontext.dof_indices);
    }

  STOP_LOG('assemble_qoi_derivative()', 'FEMSystem');
}
 

void FEMSystem::assembly (boolget_residual, boolget_jacobian) [virtual, inherited]Prepares matrix or rhs for matrix assembly. Users may reimplement this to add pre- or post-assembly code before or after calling FEMSystem::assembly()

Implements DifferentiableSystem.

Definition at line 78 of file fem_system.C.

References MeshBase::active_local_elements_begin(), MeshBase::active_local_elements_end(), DenseMatrix< T >::add(), SparseMatrix< T >::add_matrix(), NumericVector< T >::add_vector(), FEMSystem::build_context(), SparseMatrix< T >::close(), NumericVector< T >::close(), DifferentiableSystem::compute_internal_sides, DofMap::constrain_element_matrix(), DofMap::constrain_element_matrix_and_vector(), DofMap::constrain_element_vector(), DiffContext::dof_indices, FEMContext::elem, DiffContext::elem_jacobian, DiffContext::elem_residual, FEMContext::elem_side_fe_reinit(), AutoPtr< Tp >::get(), System::get_dof_map(), System::get_mesh(), DofObject::id(), SparseMatrix< T >::l1_norm(), NumericVector< T >::l1_norm(), DenseMatrix< T >::l1_norm(), ImplicitSystem::matrix, std::max(), Elem::n_sides(), Elem::neighbor(), FEMSystem::numerical_elem_jacobian(), FEMSystem::numerical_side_jacobian(), DifferentiableSystem::print_element_jacobians, DifferentiableSystem::print_jacobian_norms, DifferentiableSystem::print_jacobians, DifferentiableSystem::print_residual_norms, DifferentiableSystem::print_residuals, DifferentiableSystem::print_solution_norms, DifferentiableSystem::print_solutions, FEMContext::reinit(), ExplicitSystem::rhs, FEMContext::side, System::solution, DenseMatrix< T >::swap(), DifferentiableSystem::time_solver, System::update(), FEMSystem::verify_analytic_jacobians, DenseMatrix< T >::zero(), NumericVector< T >::zero(), and SparseMatrix< T >::zero().

Referenced by continuation_solve(), and solve_tangent().

{
  libmesh_assert(get_residual || get_jacobian);
  std::string log_name;
  if (get_residual && get_jacobian)
    log_name = 'assembly()';
  else if (get_residual)
    log_name = 'assembly(get_residual)';
  else
    log_name = 'assembly(get_jacobian)';

  START_LOG(log_name, 'FEMSystem');

  const MeshBase& mesh = this->get_mesh();

//  this->get_vector('_nonlinear_solution').localize
//    (*current_local_nonlinear_solution,
//     dof_map.get_send_list());
  this->update();

  if (print_solution_norms)
    {
//      this->get_vector('_nonlinear_solution').close();
      this->solution->close();
      std::cout << '|U| = '
//                << this->get_vector('_nonlinear_solution').l1_norm()
                << this->solution->l1_norm()
                << std::endl;
    }
  if (print_solutions)
    {
      unsigned int old_precision = std::cout.precision();
      std::cout.precision(16);
//      std::cout << 'U = [' << this->get_vector('_nonlinear_solution')
      std::cout << 'U = [' << *(this->solution)
                << '];' << std::endl;
      std::cout.precision(old_precision);
    }

  // Is this definitely necessary? [RHS]
  if (get_jacobian)
    matrix->zero();
  if (get_residual)
    rhs->zero();

  // Stupid C++ lets you set *Real* verify_analytic_jacobians = true!
  if (verify_analytic_jacobians > 0.5)
    {
      std::cerr << 'WARNING!  verify_analytic_jacobians was set '
                << 'to absurdly large value of '
                << verify_analytic_jacobians << std::endl;
      std::cerr << 'Resetting to 1e-6!' << std::endl;
      verify_analytic_jacobians = 1e-6;
    }

  // In time-dependent problems, the nonlinear function we're trying
  // to solve at each timestep may depend on the particular solver
  // we're using
  libmesh_assert (time_solver.get() != NULL);

  AutoPtr<DiffContext> con = this->build_context();
  FEMContext &_femcontext = libmesh_cast_ref<FEMContext&>(*con);

  // Build the residual and jacobian contributions on every active
  // mesh element on this processor
  MeshBase::const_element_iterator el =
    mesh.active_local_elements_begin();
  const MeshBase::const_element_iterator end_el =
    mesh.active_local_elements_end();

  for ( ; el != end_el; ++el)
    {
      _femcontext.reinit(*this, *el);

      bool jacobian_computed =
        time_solver->element_residual(get_jacobian, _femcontext);

      // Compute a numeric jacobian if we have to
      if (get_jacobian && !jacobian_computed)
        {
          // Make sure we didn't compute a jacobian and lie about it
          libmesh_assert(_femcontext.elem_jacobian.l1_norm() == 0.0);
          // Logging of numerical jacobians is done separately
          this->numerical_elem_jacobian(_femcontext);
        }

      // Compute a numeric jacobian if we're asked to verify the
      // analytic jacobian we got
      if (get_jacobian && jacobian_computed &&
          verify_analytic_jacobians != 0.0)
        {
          DenseMatrix<Number> analytic_jacobian(_femcontext.elem_jacobian);

          _femcontext.elem_jacobian.zero();
          // Logging of numerical jacobians is done separately
          this->numerical_elem_jacobian(_femcontext);

          Real analytic_norm = analytic_jacobian.l1_norm();
          Real numerical_norm = _femcontext.elem_jacobian.l1_norm();

          // If we can continue, we'll probably prefer the analytic jacobian
          analytic_jacobian.swap(_femcontext.elem_jacobian);

          // The matrix 'analytic_jacobian' will now hold the error matrix
          analytic_jacobian.add(-1.0, _femcontext.elem_jacobian);
          Real error_norm = analytic_jacobian.l1_norm();

          Real relative_error = error_norm /
                                std::max(analytic_norm, numerical_norm);

          if (relative_error > verify_analytic_jacobians)
            {
              std::cerr << 'Relative error ' << relative_error
                        << ' detected in analytic jacobian on element '
                        << _femcontext.elem->id() << '!' << std::endl;

              unsigned int old_precision = std::cout.precision();
              std::cout.precision(16);
              std::cout << 'J_analytic ' << _femcontext.elem->id() << ' = '
                        << _femcontext.elem_jacobian << std::endl;
              analytic_jacobian.add(1.0, _femcontext.elem_jacobian);
              std::cout << 'J_numeric ' << _femcontext.elem->id() << ' = '
                        << analytic_jacobian << std::endl;

              std::cout.precision(old_precision);

              libmesh_error();
            }
        }

      for (_femcontext.side = 0;
           _femcontext.side != _femcontext.elem->n_sides();
           ++_femcontext.side)
        {
          // Don't compute on non-boundary sides unless requested
          if (!compute_internal_sides && 
              _femcontext.elem->neighbor(_femcontext.side) != NULL)
            continue;

          // Any mesh movement has already been done (and restored,
          // if the TimeSolver isn't broken), but
          // reinitializing the side FE objects is still necessary
          _femcontext.elem_side_fe_reinit();

          DenseMatrix<Number> old_jacobian;
          // If we're in DEBUG mode, we should always verify that the
          // user's side_residual function doesn't alter our existing
          // jacobian and then lie about it
#ifndef DEBUG
          // Even if we're not in DEBUG mode, when we're verifying
          // analytic jacobians we'll want to verify each side's
          // jacobian contribution separately
          if (verify_analytic_jacobians != 0.0 && get_jacobian)
#endif // ifndef DEBUG
            {
              old_jacobian = _femcontext.elem_jacobian;
              _femcontext.elem_jacobian.zero();
            }
          jacobian_computed =
            time_solver->side_residual(get_jacobian, _femcontext);

          // Compute a numeric jacobian if we have to
          if (get_jacobian && !jacobian_computed)
            {
              // In DEBUG mode, we've already set elem_jacobian == 0,
              // so we can make sure side_residual didn't compute a
              // jacobian and lie about it
#ifdef DEBUG
              libmesh_assert(_femcontext.elem_jacobian.l1_norm() == 0.0);
#endif
              // Logging of numerical jacobians is done separately
              this->numerical_side_jacobian(_femcontext);

              // If we're in DEBUG mode or if
              // verify_analytic_jacobians is on, we've moved
              // elem_jacobian's accumulated values into old_jacobian.
              // Now let's add them back.
#ifndef DEBUG
              if (verify_analytic_jacobians != 0.0)
#endif // ifndef DEBUG
                _femcontext.elem_jacobian += old_jacobian;
            }

          // Compute a numeric jacobian if we're asked to verify the
          // analytic jacobian we got
          if (get_jacobian && jacobian_computed &&
              verify_analytic_jacobians != 0.0)
            {
              DenseMatrix<Number> analytic_jacobian(_femcontext.elem_jacobian);

              _femcontext.elem_jacobian.zero();
              // Logging of numerical jacobians is done separately
              this->numerical_side_jacobian(_femcontext);

              Real analytic_norm = analytic_jacobian.l1_norm();
              Real numerical_norm = _femcontext.elem_jacobian.l1_norm();

              // If we can continue, we'll probably prefer the analytic jacobian
              analytic_jacobian.swap(_femcontext.elem_jacobian);

              // The matrix 'analytic_jacobian' will now hold the error matrix
              analytic_jacobian.add(-1.0, _femcontext.elem_jacobian);
              Real error_norm = analytic_jacobian.l1_norm();

              Real relative_error = error_norm /
                                    std::max(analytic_norm, numerical_norm);

              if (relative_error > verify_analytic_jacobians)
                {
                  std::cerr << 'Relative error ' << relative_error
                            << ' detected in analytic jacobian on element '
                            << _femcontext.elem->id()
                            << ', side '
                            << _femcontext.side << '!' << std::endl;

                  unsigned int old_precision = std::cout.precision();
                  std::cout.precision(16);
                  std::cout << 'J_analytic ' << _femcontext.elem->id() << ' = '
                            << _femcontext.elem_jacobian << std::endl;
                  analytic_jacobian.add(1.0, _femcontext.elem_jacobian);
                  std::cout << 'J_numeric ' << _femcontext.elem->id() << ' = '
                            << analytic_jacobian << std::endl;
                  std::cout.precision(old_precision);

                  libmesh_error();
                }
              // Once we've verified a side, we'll want to add back the
              // rest of the accumulated jacobian
              _femcontext.elem_jacobian += old_jacobian;
            }
          // In DEBUG mode, we've set elem_jacobian == 0, and we
          // may still need to add the old jacobian back
#ifdef DEBUG
          if (get_jacobian && jacobian_computed &&
              verify_analytic_jacobians == 0.0)
            {
              _femcontext.elem_jacobian += old_jacobian;
            }
#endif // ifdef DEBUG
        }

#ifdef LIBMESH_ENABLE_AMR
      // We turn off the asymmetric constraint application;
      // enforce_constraints_exactly() should be called in the solver
      if (get_residual && get_jacobian)
        this->get_dof_map().constrain_element_matrix_and_vector
          (_femcontext.elem_jacobian, _femcontext.elem_residual,
           _femcontext.dof_indices, false);
      else if (get_residual)
        this->get_dof_map().constrain_element_vector
          (_femcontext.elem_residual, _femcontext.dof_indices, false);
      else if (get_jacobian)
        this->get_dof_map().constrain_element_matrix
          (_femcontext.elem_jacobian, _femcontext.dof_indices, false);
#endif // #ifdef LIBMESH_ENABLE_AMR

      if (get_jacobian && print_element_jacobians)
        {
          unsigned int old_precision = std::cout.precision();
          std::cout.precision(16);
          std::cout << 'J_elem ' << _femcontext.elem->id() << ' = '
                    << _femcontext.elem_jacobian << std::endl;
          std::cout.precision(old_precision);
        }

      if (get_jacobian)
        this->matrix->add_matrix (_femcontext.elem_jacobian,
                                  _femcontext.dof_indices);
      if (get_residual)
        this->rhs->add_vector (_femcontext.elem_residual,
                               _femcontext.dof_indices);
    }


  if (get_residual && print_residual_norms)
    {
      this->rhs->close();
      std::cout << '|F| = ' << this->rhs->l1_norm() << std::endl;
    }
  if (get_residual && print_residuals)
    {
      this->rhs->close();
      unsigned int old_precision = std::cout.precision();
      std::cout.precision(16);
      std::cout << 'F = [' << *(this->rhs) << '];' << std::endl;
      std::cout.precision(old_precision);
    }
  if (get_jacobian && print_jacobian_norms)
    {
      this->matrix->close();
      std::cout << '|J| = ' << this->matrix->l1_norm() << std::endl;
    }
  if (get_jacobian && print_jacobians)
    {
      this->matrix->close();
      unsigned int old_precision = std::cout.precision();
      std::cout.precision(16);
      std::cout << 'J = [' << *(this->matrix) << '];' << std::endl;
      std::cout.precision(old_precision);
    }
  STOP_LOG(log_name, 'FEMSystem');
}
 

void System::attach_assemble_function (void fptrEquationSystems &es,const std::string &name) [inherited]Register a user function to use in assembling the system matrix and RHS.

Definition at line 1084 of file system.C.

References System::_assemble_system.

{
  libmesh_assert (fptr != NULL);
  
  _assemble_system = fptr;  
}
 

void System::attach_constraint_function (void fptrEquationSystems &es,const std::string &name) [inherited]Register a user function for imposing constraints.

Definition at line 1094 of file system.C.

References System::_constrain_system.

{
  libmesh_assert (fptr != NULL);
  
  _constrain_system = fptr;  
}
 

void System::attach_init_function (void fptrEquationSystems &es,const std::string &name) [inherited]Register a user function to use in initializing the system.

Definition at line 1074 of file system.C.

References System::_init_system.

{
  libmesh_assert (fptr != NULL);
  
  _init_system = fptr;
}
 

void System::attach_QOI_derivative (void fptrEquationSystems &es,const std::string &name) [inherited]Register a user function for evaluating derivatives of a quantity of interest with respect to test functions, whose values should be placed in System::rhs

Definition at line 1114 of file system.C.

References System::_qoi_evaluate_derivative.

{
  libmesh_assert (fptr != NULL);
  
  _qoi_evaluate_derivative = fptr;  
}
 

void System::attach_QOI_function (void fptrEquationSystems &es,const std::string &name) [inherited]Register a user function for evaluating a quantity of interest, whose value should be placed in System::qoi

Definition at line 1104 of file system.C.

References System::_qoi_evaluate.

{
  libmesh_assert (fptr != NULL);
  
  _qoi_evaluate = fptr;  
}
 

AutoPtr< DiffContext > FEMSystem::build_context () [virtual, inherited]Builds a FEMContext object with enough information to do evaluations on each element.

For most problems, the default FEMSystem implementation is correct; users who subclass FEMContext will need to also reimplement this method to build it.

Reimplemented from DifferentiableSystem.

Definition at line 724 of file fem_system.C.

Referenced by FEMSystem::assemble_qoi(), FEMSystem::assemble_qoi_derivative(), FEMSystem::assembly(), FEMSystem::mesh_position_get(), FEMSystem::mesh_position_set(), and FEMSystem::postprocess().

{
  return AutoPtr<DiffContext>(new FEMContext(*this));
}
 

Real System::calculate_norm (NumericVector< Number > &v, unsigned intvar = 0, FEMNormTypenorm_type = L2) const [inherited]Returns:

a norm of variable var in the vector v, in the specified norm (e.g. L2, H0, H1)

Definition at line 824 of file system.C.

References libMeshEnums::DISCRETE_L1, libMeshEnums::DISCRETE_L2, libMeshEnums::DISCRETE_L_INF, System::discrete_var_norm(), libMeshEnums::L2, and System::n_vars().

Referenced by AdaptiveTimeSolver::calculate_norm(), and UnsteadySolver::du().

{
  //short circuit to save time
  if(norm_type == DISCRETE_L1 ||
     norm_type == DISCRETE_L2 ||
     norm_type == DISCRETE_L_INF)
    return discrete_var_norm(v,var,norm_type);

  // Not a discrete norm
  std::vector<FEMNormType> norms(this->n_vars(), L2);
  std::vector<Real> weights(this->n_vars(), 0.0);
  norms[var] = norm_type;
  weights[var] = 1.0;
  Real val = this->calculate_norm(v, SystemNorm(norms, weights));
  return val;
}
 

Real System::calculate_norm (NumericVector< Number > &v, const SystemNorm &norm) const [inherited]Returns:

a norm of the vector v, using component_norm and component_scale to choose and weight the norms of each variable.

Definition at line 845 of file system.C.

References System::_dof_map, MeshBase::active_local_elements_begin(), MeshBase::active_local_elements_end(), TypeTensor< T >::add_scaled(), TypeVector< T >::add_scaled(), FEBase::build(), FEType::default_quadrature_rule(), libMeshEnums::DISCRETE_L1, libMeshEnums::DISCRETE_L2, libMeshEnums::DISCRETE_L_INF, System::discrete_var_norm(), DofMap::dof_indices(), AutoPtr< Tp >::get(), System::get_dof_map(), System::get_mesh(), libMeshEnums::H1, libMeshEnums::H1_SEMINORM, libMeshEnums::H2, libMeshEnums::H2_SEMINORM, SystemNorm::is_discrete(), NumericVector< T >::l1_norm(), libMeshEnums::L2, NumericVector< T >::l2_norm(), libmesh_norm(), NumericVector< T >::linfty_norm(), NumericVector< T >::localize(), MeshBase::mesh_dimension(), System::n_vars(), libMeshEnums::SERIAL, NumericVector< T >::size(), TypeTensor< T >::size_sq(), TypeVector< T >::size_sq(), SystemNorm::type(), DofMap::variable_type(), SystemNorm::weight(), and SystemNorm::weight_sq().

{
  // This function must be run on all processors at once
  parallel_only();

  START_LOG ('calculate_norm()', 'System');

  // Zero the norm before summation
  Real v_norm = 0.;

  if (norm.is_discrete())
    {
      STOP_LOG ('calculate_norm()', 'System');
      //Check to see if all weights are 1.0
      unsigned int check_var = 0;
      for (; check_var != this->n_vars(); ++check_var)
        if(norm.weight(check_var) != 1.0)
          break;

      //All weights were 1.0 so just do the full vector discrete norm
      if(check_var == this->n_vars())
        {
          FEMNormType norm_type = norm.type(0);
          
          if(norm_type == DISCRETE_L1)
            return v.l1_norm();
          if(norm_type == DISCRETE_L2)
            return v.l2_norm();
          if(norm_type == DISCRETE_L_INF)
            return v.linfty_norm();
          else
            libmesh_error();
        }

      for (unsigned int var=0; var != this->n_vars(); ++var)
        {
          // Skip any variables we don't need to integrate
          if (norm.weight(var) == 0.0)
            continue;

          v_norm += norm.weight(var) * discrete_var_norm(v, var, norm.type(var));
        }

      return v_norm;
    }

  // Localize the potentially parallel vector
  AutoPtr<NumericVector<Number> > local_v = NumericVector<Number>::build();
  local_v->init(v.size(), true, SERIAL);
  v.localize (*local_v, _dof_map->get_send_list()); 

  unsigned int dim = this->get_mesh().mesh_dimension();

  // Loop over all variables
  for (unsigned int var=0; var != this->n_vars(); ++var)
    {
      // Skip any variables we don't need to integrate
      if (norm.weight(var) == 0.0)
        continue;

      const FEType& fe_type = this->get_dof_map().variable_type(var);
      AutoPtr<QBase> qrule =
        fe_type.default_quadrature_rule (dim);
      AutoPtr<FEBase> fe
        (FEBase::build(dim, fe_type));
      fe->attach_quadrature_rule (qrule.get());

      const std::vector<Real>&               JxW = fe->get_JxW();
      const std::vector<std::vector<Real> >* phi = NULL;
      if (norm.type(var) == H1 ||
          norm.type(var) == H2 ||
          norm.type(var) == L2)
        phi = &(fe->get_phi());

      const std::vector<std::vector<RealGradient> >* dphi = NULL;
      if (norm.type(var) == H1 ||
          norm.type(var) == H2 ||
          norm.type(var) == H1_SEMINORM)
        dphi = &(fe->get_dphi());
#ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
      const std::vector<std::vector<RealTensor> >*   d2phi = NULL;
      if (norm.type(var) == H2 ||
          norm.type(var) == H2_SEMINORM)
        d2phi = &(fe->get_d2phi());
#endif

      std::vector<unsigned int> dof_indices;

      // Begin the loop over the elements
      MeshBase::const_element_iterator       el     =
        this->get_mesh().active_local_elements_begin();
      const MeshBase::const_element_iterator end_el =
        this->get_mesh().active_local_elements_end();

      for ( ; el != end_el; ++el)
        {
          const Elem* elem = *el;

          fe->reinit (elem);

          this->get_dof_map().dof_indices (elem, dof_indices, var);

          const unsigned int n_qp = qrule->n_points();

          const unsigned int n_sf = dof_indices.size();

          // Begin the loop over the Quadrature points.
          for (unsigned int qp=0; qp<n_qp; qp++)
            {
              if (norm.type(var) == H1 ||
                  norm.type(var) == H2 ||
                  norm.type(var) == L2)
                {
                  Number u_h = 0.;
                  for (unsigned int i=0; i != n_sf; ++i)
                    u_h += (*phi)[i][qp] * (*local_v)(dof_indices[i]);
                  v_norm += norm.weight_sq(var) *
                            JxW[qp] * libmesh_norm(u_h);
                }

              if (norm.type(var) == H1 ||
                  norm.type(var) == H2 ||
                  norm.type(var) == H1_SEMINORM)
                {
                  Gradient grad_u_h;
                  for (unsigned int i=0; i != n_sf; ++i)
                    grad_u_h.add_scaled((*dphi)[i][qp], (*local_v)(dof_indices[i]));
                  v_norm += norm.weight_sq(var) *
                            JxW[qp] * grad_u_h.size_sq();
                }

#ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
              if (norm.type(var) == H2 ||
                  norm.type(var) == H2_SEMINORM)
                {
                  Tensor hess_u_h;
                  for (unsigned int i=0; i != n_sf; ++i)
                    hess_u_h.add_scaled((*d2phi)[i][qp], (*local_v)(dof_indices[i]));
                  v_norm += norm.weight_sq(var) *
                            JxW[qp] * hess_u_h.size_sq();
                }
#endif
            }
        }
    }

  Parallel::sum(v_norm);

  STOP_LOG ('calculate_norm()', 'System');

  return std::sqrt(v_norm);
}
 

void ContinuationSystem::clear () [virtual]Clear all the data structures associated with the system.

Reimplemented from FEMSystem.

Definition at line 74 of file continuation_system.C.

References FEMSystem::clear().

Referenced by ~ContinuationSystem().

{
  // FIXME: Do anything here, e.g. zero vectors, etc?

  // Call the Parent's clear function
  Parent::clear();
}
 

bool System::compare (const System &other_system, const Realthreshold, const boolverbose) const [virtual, inherited]Returns:

true when the other system contains identical data, up to the given threshold. Outputs some diagnostic info when verbose is set.

Definition at line 380 of file system.C.

References System::_can_add_vectors, System::_sys_name, System::_vectors, AutoPtr< Tp >::get(), System::get_vector(), System::n_vectors(), System::name(), and System::solution.

Referenced by EquationSystems::compare().

{
  // we do not care for matrices, but for vectors
  libmesh_assert (!_can_add_vectors);
  libmesh_assert (!other_system._can_add_vectors);

  if (verbose)
    {
      std::cout << '  Systems '' << _sys_name << ''' << std::endl;
      std::cout << '   comparing matrices not supported.' << std::endl;
      std::cout << '   comparing names...';
    }

  // compare the name: 0 means identical
  const int name_result = _sys_name.compare(other_system.name());
  if (verbose)
    {
      if (name_result == 0)
        std::cout << ' identical.' << std::endl;
      else
        std::cout << '  names not identical.' << std::endl;
      std::cout << '   comparing solution vector...';
    }


  // compare the solution: -1 means identical
  const int solu_result = solution->compare (*other_system.solution.get(),
                                             threshold);

  if (verbose)
    {
      if (solu_result == -1)
        std::cout << ' identical up to threshold.' << std::endl;
      else
        std::cout << '  first difference occured at index = ' 
                  << solu_result << '.' << std::endl;
    }


  // safety check, whether we handle at least the same number
  // of vectors
  std::vector<int> ov_result;

  if (this->n_vectors() != other_system.n_vectors())
    {
      if (verbose)
        {
          std::cout << '   Fatal difference. This system handles ' 
                    << this->n_vectors() << ' add'l vectors,' << std::endl
                    << '   while the other system handles '
                    << other_system.n_vectors() 
                    << ' add'l vectors.' << std::endl
                    << '   Aborting comparison.' << std::endl;
        }
      return false;
    }
  else if (this->n_vectors() == 0)
    {
      // there are no additional vectors...
      ov_result.clear ();
    }
  else
    {
      // compare other vectors
      for (const_vectors_iterator pos = _vectors.begin();
           pos != _vectors.end(); ++pos)
        {
          if (verbose)
              std::cout << '   comparing vector ''
                        << pos->first << '' ...';

          // assume they have the same name
          const NumericVector<Number>& other_system_vector = 
              other_system.get_vector(pos->first);

          ov_result.push_back(pos->second->compare (other_system_vector,
                                                    threshold));

          if (verbose)
            {
              if (ov_result[ov_result.size()-1] == -1)
                std::cout << ' identical up to threshold.' << std::endl;
              else
                std::cout << ' first difference occured at' << std::endl
                          << '   index = ' << ov_result[ov_result.size()-1] << '.' << std::endl;
            }

        }

    } // finished comparing additional vectors


  bool overall_result;
       
  // sum up the results
  if ((name_result==0) && (solu_result==-1))
    {
      if (ov_result.size()==0)
        overall_result = true;
      else
        {
          bool ov_identical;
          unsigned int n    = 0;
          do
            {
              ov_identical = (ov_result[n]==-1);
              n++;
            }
          while (ov_identical && n<ov_result.size());
          overall_result = ov_identical;
        }
    }
  else
    overall_result = false;

  if (verbose)
    {
      std::cout << '   finished comparisons, ';
      if (overall_result)
        std::cout << 'found no differences.' << std::endl << std::endl;
      else 
        std::cout << 'found differences.' << std::endl << std::endl;
    }
          
  return overall_result;
}
 

void ContinuationSystem::continuation_solve ()Perform a continuation solve of the system. In general, you can only begin the continuation solves after either reading in or solving for two previous values of the control parameter. The prior two solutions are required for starting up the continuation method.

Definition at line 352 of file continuation_system.C.

References NumericVector< T >::add(), apply_predictor(), FEMSystem::assembly(), SparseMatrix< T >::close(), NumericVector< T >::close(), continuation_parameter, continuation_parameter_tolerance, delta_u, dlambda_ds, NumericVector< T >::dot(), ds_current, du_ds, G_Lambda, AutoPtr< Tp >::get(), initial_newton_tolerance, initialize_tangent(), NumericVector< T >::l2_norm(), libmesh_real(), linear_solver, ImplicitSystem::matrix, max_continuation_parameter, DiffSolver::max_linear_iterations, DiffSolver::max_nonlinear_iterations, std::min(), min_continuation_parameter, n_arclength_reductions, n_backtrack_steps, newton_progress_check, newton_solver, newton_step, old_continuation_parameter, Utility::pow(), previous_u, quiet, Residual, ExplicitSystem::rhs, rhs_mode, NumericVector< T >::scale(), System::solution, solution_tolerance, tangent_initialized, Theta, Theta_LOCA, DifferentiableSystem::time_solver, y, y_old, z, and NumericVector< T >::zero().

{
  // Be sure the user has set the continuation parameter pointer
  if (!continuation_parameter)
    {
      std::cerr << 'You must set the continuation_parameter pointer '
                << 'to a member variable of the derived class, preferably in the '
                << 'Derived class's init_data function.  This is how the ContinuationSystem '
                << 'updates the continuation parameter.'
                << std::endl;
      
      libmesh_error();
    }

  // Use extra precision for all the numbers printed in this function.
  unsigned int old_precision = std::cout.precision();
  std::cout.precision(16);
  std::cout.setf(std::ios_base::scientific);
  
  // We can't start solving the augmented PDE system unless the tangent
  // vectors have been initialized.  This only needs to occur once.
  if (!tangent_initialized)
    initialize_tangent();

  // Save the old value of -du/dlambda.  This will be used after the Newton iterations
  // to compute the angle between previous tangent vectors.  This cosine of this angle is
  //
  // tau := abs( (du/d(lambda)_i , du/d(lambda)_{i-1}) / (||du/d(lambda)_i|| * ||du/d(lambda)_{i-1}||) )
  //
  // The scaling factor tau (which should vary between 0 and 1) is used to shrink the step-size ds
  // when we are approaching a turning point.  Note that it can only shrink the step size.  
  *y_old = *y;
  
  // Set pointer to underlying Newton solver
  if (!newton_solver)
    newton_solver = libmesh_cast_ptr<NewtonSolver*> (this->time_solver->diff_solver().get());
  
  // A pair for catching return values from linear system solves.
  std::pair<unsigned int, Real> rval;

  // Convergence flag for the entire arcstep
  bool arcstep_converged = false;

  // Begin loop over arcstep reductions.
  for (unsigned int ns=0; ns<n_arclength_reductions; ++ns)
    {
      if (!quiet)
        {
          std::cout << 'Current arclength stepsize, ds_current=' << ds_current << std::endl;
          std::cout << 'Current parameter value, lambda=' << *continuation_parameter << std::endl;
        }
      
      // Upon exit from the nonlinear loop, the newton_converged flag
      // will tell us the convergence status of Newton's method.
      bool newton_converged = false;

      // The nonlinear residual before *any* nonlinear steps have been taken.
      Real nonlinear_residual_firststep = 0.;

      // The nonlinear residual from the current 'k' Newton step, before the Newton step
      Real nonlinear_residual_beforestep = 0.;

      // The nonlinear residual from the current 'k' Newton step, after the Newton step
      Real nonlinear_residual_afterstep = 0.;

      // The linear solver tolerance, can be updated dynamically at each Newton step.
      Real current_linear_tolerance = 0.;
      
      // The nonlinear loop
      for (newton_step=0; newton_step<newton_solver->max_nonlinear_iterations; ++newton_step)
        {
          std::cout << '=== Starting Newton step ' << newton_step << ' ===' << std::endl;
          
          // Set the linear system solver tolerance
//        // 1.) Set the current linear tolerance based as a multiple of the current residual of the system.
//        const Real residual_multiple = 1.e-4;
//        Real current_linear_tolerance = residual_multiple*nonlinear_residual_beforestep;

//        // But if the current residual isn't small, don't let the solver exit with zero iterations!
//        if (current_linear_tolerance > 1.)
//          current_linear_tolerance = residual_multiple;

          // 2.) Set the current linear tolerance based on the method based on technique of Eisenstat & Walker.
          if (newton_step==0)
            {
              // At first step, only try reducing the residual by a small amount
              current_linear_tolerance = initial_newton_tolerance;//0.01;
            }

          else
            {
              // The new tolerance is based on the ratio of the most recent tolerances
              const Real alp=0.5*(1.+std::sqrt(5.));
              const Real gam=0.9;

              libmesh_assert (nonlinear_residual_beforestep != 0.0);
              libmesh_assert (nonlinear_residual_afterstep != 0.0);

              current_linear_tolerance = std::min(gam*std::pow(nonlinear_residual_afterstep/nonlinear_residual_beforestep, alp),
                                                  current_linear_tolerance*current_linear_tolerance
                                                  );

              // Don't let it get ridiculously small!!
              if (current_linear_tolerance < 1.e-12)
                current_linear_tolerance = 1.e-12;
            }

          if (!quiet)
            std::cout << 'Using current_linear_tolerance=' << current_linear_tolerance << std::endl;

          
          // Assemble the residual (and Jacobian).
          rhs_mode = Residual;
          assembly(true,   // Residual
                   true); // Jacobian
          rhs->close();

          // Save the current nonlinear residual.  We don't need to recompute the residual unless
          // this is the first step, since it was already computed as part of the convergence check
          // at the end of the last loop iteration.
          if (newton_step==0)
            {
              nonlinear_residual_beforestep = rhs->l2_norm();

              // Store the residual before any steps have been taken.  This will *not*
              // be updated at each step, and can be used to see if any progress has
              // been made from the initial residual at later steps.
              nonlinear_residual_firststep = nonlinear_residual_beforestep;

              const Real old_norm_u = solution->l2_norm();
              std::cout << '  (before step) ||R||_{L2} = ' << nonlinear_residual_beforestep << std::endl;
              std::cout << '  (before step) ||R||_{L2}/||u|| = ' << nonlinear_residual_beforestep / old_norm_u << std::endl;

              // In rare cases (very small arcsteps), it's possible that the residual is
              // already below our absolute linear tolerance.
              if (nonlinear_residual_beforestep  < solution_tolerance)
                {
                  if (!quiet)
                    std::cout << 'Initial guess satisfied linear tolerance, exiting with zero Newton iterations!' << std::endl;

                  // Since we go straight from here to the solve of the next tangent, we
                  // have to close the matrix before it can be assembled again.
                  matrix->close();
                  newton_converged=true;
                  break; // out of Newton iterations, with newton_converged=true
                }
            }

          else
            {
              nonlinear_residual_beforestep = nonlinear_residual_afterstep;
            }

          
          // Solve the linear system G_u*z = G
          // Initial guess?
          z->zero(); // It seems to be extremely important to zero z here, otherwise the solver quits early.
          z->close();
          
          // It's possible that we have selected the current_linear_tolerance so large that
          // a guess of z=zero yields a linear system residual |Az + R| small enough that the
          // linear solver exits in zero iterations.  If this happens, we will reduce the
          // current_linear_tolerance until the linear solver does at least 1 iteration.
          do
            {
              rval =
                linear_solver->solve(*matrix,
                                     *z,
                                     *rhs,
                                     //1.e-12,
                                     current_linear_tolerance,
                                     newton_solver->max_linear_iterations);   // max linear iterations
              
              if (rval.first==0)
                {
                  if (newton_step==0)
                    {
                      std::cout << 'Repeating initial solve with smaller linear tolerance!' << std::endl;
                      current_linear_tolerance *= initial_newton_tolerance; // reduce the linear tolerance to force the solver to do some work
                    }
                  else
                    {
                      // We shouldn't get here ... it means the linear solver did no work on a Newton
                      // step other than the first one.  If this happens, we need to think more about our
                      // tolerance selection.
                      libmesh_error();
                    }
                }
              
            } while (rval.first==0);

          
          if (!quiet)
            std::cout << '  G_u*z = G solver converged at step '
                      << rval.first
                      << ' linear tolerance = '
                      << rval.second
                      << '.'
                      << std::endl;

          // Sometimes (I am not sure why) the linear solver exits after zero iterations.
          // Perhaps it is hitting PETSc's divergence tolerance dtol???  If this occurs,
          // we should break out of the Newton iteration loop because nothing further is
          // going to happen...  Of course if the tolerance is already small enough after
          // zero iterations (how can this happen?!) we should not quit.
          if ((rval.first == 0) && (rval.second > current_linear_tolerance*nonlinear_residual_beforestep))
            {
              if (!quiet)
                std::cout << 'Linear solver exited in zero iterations!' << std::endl;

              // Try to find out the reason for convergence/divergence
              linear_solver->print_converged_reason();
              
              break; // out of Newton iterations
            }
          
          // Note: need to scale z by -1 since our code always solves Jx=R
          // instead of Jx=-R.
          z->scale(-1.);
          z->close();





          
          // Assemble the G_Lambda vector, skip residual.
          rhs_mode = G_Lambda;

          // Assemble both rhs and Jacobian
          assembly(true,  // Residual
                   false); // Jacobian

          // Not sure if this is really necessary
          rhs->close();
          const Real yrhsnorm=rhs->l2_norm();
          if (yrhsnorm == 0.0)
            {
              std::cout << '||G_Lambda|| = 0' << std::endl;
              libmesh_error();
            }

          // We select a tolerance for the y-system which is based on the inexact Newton
          // tolerance but scaled by an extra term proportional to the RHS (which is not -> 0 in this case)
          const Real ysystemtol=current_linear_tolerance*(nonlinear_residual_beforestep/yrhsnorm);
          if (!quiet)
            std::cout << 'ysystemtol=' << ysystemtol << std::endl;
          
          // Solve G_u*y = G_{
} // FIXME: Initial guess? This is really a solve for -du/dlambda so we could try // initializing it with the latest approximation to that... du/dlambda ~ du/ds * ds/dlambda //*y = *solution; //y->add(-1., *previous_u); //y->scale(-1. / (*continuation_parameter - old_continuation_parameter)); // Be careful of divide by zero... //y->close(); // const unsigned int max_attempts=1; // unsigned int attempt=0; // do // { // if (!quiet) // std::cout << 'Trying to solve tangent system, attempt ' << attempt << std::endl; rval = linear_solver->solve(*matrix, *y, *rhs, //1.e-12, ysystemtol, newton_solver->max_linear_iterations); // max linear iterations if (!quiet) std::cout << ' G_u*y = G_{lambda} solver converged at step ' << rval.first << ', linear tolerance = ' << rval.second << '.' << std::endl; // Sometimes (I am not sure why) the linear solver exits after zero iterations. // Perhaps it is hitting PETSc's divergence tolerance dtol??? If this occurs, // we should break out of the Newton iteration loop because nothing further is // going to happen... if ((rval.first == 0) && (rval.second > ysystemtol)) { if (!quiet) std::cout << 'Linear solver exited in zero iterations!' << std::endl; break; // out of Newton iterations } // ++attempt; // } while ((attempt<max_attempts) && (rval.first==newton_solver->max_linear_iterations)); // Compute N, the residual of the arclength constraint eqn. // Note 1: N(u,lambda,s) := (u-u_{old}, du_ds) + (lambda-lambda_{old}, dlambda_ds) - _ds // We temporarily use the delta_u vector as a temporary vector for this calculation. *delta_u = *solution; delta_u->add(-1., *previous_u); // First part of the arclength constraint const Number N1 = Theta_LOCA*Theta_LOCA*Theta*delta_u->dot(*du_ds); const Number N2 = ((*continuation_parameter) - old_continuation_parameter)*dlambda_ds; const Number N3 = ds_current; if (!quiet) { std::cout << ' N1=' << N1 << std::endl; std::cout << ' N2=' << N2 << std::endl; std::cout << ' N3=' << N3 << std::endl; } // The arclength constraint value const Number N = N1+N2-N3; if (!quiet) std::cout << ' N=' << N << std::endl; const Number duds_dot_z = du_ds->dot(*z); const Number duds_dot_y = du_ds->dot(*y); //std::cout << 'duds_dot_z=' << duds_dot_z << std::endl; //std::cout << 'duds_dot_y=' << duds_dot_y << std::endl; //std::cout << 'dlambda_ds=' << dlambda_ds << std::endl; const Number delta_lambda_numerator = -(N + Theta_LOCA*Theta_LOCA*Theta*duds_dot_z); const Number delta_lambda_denominator = (dlambda_ds - Theta_LOCA*Theta_LOCA*Theta*duds_dot_y); libmesh_assert (delta_lambda_denominator != 0.0); // Now, we are ready to compute the step delta_lambda const Number delta_lambda_comp = delta_lambda_numerator / delta_lambda_denominator; // Lambda is real-valued const Real delta_lambda = libmesh_real(delta_lambda_comp); // Knowing delta_lambda, we are ready to update delta_u // delta_u = z - delta_lambda*y delta_u->zero(); delta_u->add(1., *z); delta_u->add(-delta_lambda, *y); delta_u->close(); // Update the system solution and the continuation parameter. solution->add(1., *delta_u); solution->close(); *continuation_parameter += delta_lambda; // Did the Newton step actually reduce the residual? rhs_mode = Residual; assembly(true, // Residual false); // Jacobian rhs->close(); nonlinear_residual_afterstep = rhs->l2_norm(); // In a 'normal' Newton step, ||du||/||R|| > 1 since the most recent // step is where you 'just were' and the current residual is where // you are now. It can occur that ||du||/||R|| < 1, but these are // likely not good cases to attempt backtracking (?). const Real norm_du_norm_R = delta_u->l2_norm() / nonlinear_residual_afterstep; if (!quiet) std::cout << ' norm_du_norm_R=' << norm_du_norm_R << std::endl; // Factor to decrease the stepsize by for backtracking Real newton_stepfactor = 1.; const bool attempt_backtracking = (nonlinear_residual_afterstep > solution_tolerance) && (nonlinear_residual_afterstep > nonlinear_residual_beforestep) && (n_backtrack_steps>0) && (norm_du_norm_R > 1.) ; // If residual is not reduced, do Newton back tracking. if (attempt_backtracking) { if (!quiet) std::cout << 'Newton step did not reduce residual.' << std::endl; // back off the previous step. solution->add(-1., *delta_u); solution->close(); *continuation_parameter -= delta_lambda; // Backtracking: start cutting the Newton stepsize by halves until // the new residual is actually smaller... for (unsigned int backtrack_step=0; backtrack_step<n_backtrack_steps; ++backtrack_step) { newton_stepfactor *= 0.5; if (!quiet) std::cout << 'Shrinking step size by ' << newton_stepfactor << std::endl; // Take fractional step solution->add(newton_stepfactor, *delta_u); solution->close(); *continuation_parameter += newton_stepfactor*delta_lambda; rhs_mode = Residual; assembly(true, // Residual false); // Jacobian rhs->close(); nonlinear_residual_afterstep = rhs->l2_norm(); if (!quiet) std::cout << 'At shrink step ' << backtrack_step << ', nonlinear_residual_afterstep=' << nonlinear_residual_afterstep << std::endl; if (nonlinear_residual_afterstep < nonlinear_residual_beforestep) { if (!quiet) std::cout << 'Backtracking succeeded!' << std::endl; break; // out of backtracking loop } else { // Back off that step solution->add(-newton_stepfactor, *delta_u); solution->close(); *continuation_parameter -= newton_stepfactor*delta_lambda; } // Save a copy of the solution from before the Newton step. //AutoPtr<NumericVector<Number> > prior_iterate = solution->clone(); } } // end if (attempte_backtracking) // If we tried backtracking but the residual is still not reduced, print message. if ((attempt_backtracking) && (nonlinear_residual_afterstep > nonlinear_residual_beforestep)) { //std::cerr << 'Backtracking failed.' << std::endl; std::cout << 'Backtracking failed.' << std::endl; // 1.) Quit, exit program. //libmesh_error(); // 2.) Continue with last newton_stepfactor if (newton_step<3) { solution->add(newton_stepfactor, *delta_u); solution->close(); *continuation_parameter += newton_stepfactor*delta_lambda; if (!quiet) std::cout << 'Backtracking could not reduce residual ... continuing anyway!' << std::endl; } // 3.) Break out of Newton iteration loop with newton_converged = false, // reduce the arclength stepsize, and try again. else { break; // out of Newton iteration loop, with newton_converged=false } } // Another type of convergence check: suppose the residual has not been reduced // from its initial value after half of the allowed Newton steps have occurred. // In our experience, this typically means that it isn't going to converge and // we could probably save time by dropping out of the Newton iteration loop and // trying a smaller arcstep. if (this->newton_progress_check) { if ((nonlinear_residual_afterstep > nonlinear_residual_firststep) && (newton_step+1 > static_cast<unsigned int>(0.5*newton_solver->max_nonlinear_iterations))) { std::cout << 'Progress check failed: the current residual: ' << nonlinear_residual_afterstep << ', is << 'larger than the initial residual, and half of the allowed << 'number of Newton iterations have elapsed. << 'Exiting Newton iterations with converged==false.' << std::endl; break; // out of Newton iteration loop, newton_converged = false } } // Safety check: Check the current continuation parameter against user-provided min-allowable parameter value if (*continuation_parameter < min_continuation_parameter) { std::cout << 'Continuation parameter fell below min-allowable value.' << std::endl; // libmesh_error(); break; // out of Newton iteration loop, newton_converged = false } // Safety check: Check the current continuation parameter against user-provided max-allowable parameter value if ( (max_continuation_parameter != 0.0) && (*continuation_parameter > max_continuation_parameter) ) { std::cout << 'Current continuation parameter value: ' << *continuation_parameter << ' exceeded max-allowable value.' << std::endl; // libmesh_error(); break; // out of Newton iteration loop, newton_converged = false } // Check the convergence of the parameter and the solution. If they are small // enough, we can break out of the Newton iteration loop. const Real norm_delta_u = delta_u->l2_norm(); const Real norm_u = solution->l2_norm(); std::cout << ' delta_lambda = ' << delta_lambda << std::endl; std::cout << ' newton_stepfactor*delta_lambda = ' << newton_stepfactor*delta_lambda << std::endl; std::cout << ' lambda_current = ' << *continuation_parameter << std::endl; std::cout << ' ||delta_u|| = ' << norm_delta_u << std::endl; std::cout << ' ||delta_u||/||u|| = ' << norm_delta_u / norm_u << std::endl; // Evaluate the residual at the current Newton iterate. We don't want to detect // convergence due to a small Newton step when the residual is still not small. rhs_mode = Residual; assembly(true, // Residual false); // Jacobian rhs->close(); const Real norm_residual = rhs->l2_norm(); std::cout << ' ||R||_{L2} = ' << norm_residual << std::endl; std::cout << ' ||R||_{L2}/||u|| = ' << norm_residual / norm_u << std::endl; // FIXME: The norm_delta_u tolerance (at least) should be relative. // It doesn't make sense to converge a solution whose size is ~ 10^5 to // a tolerance of 1.e-6. Oh, and we should also probably check the // (relative) size of the residual as well, instead of just the step. if ((std::abs(delta_lambda) < continuation_parameter_tolerance) && //(norm_delta_u < solution_tolerance) && // This is a *very* strict criterion we can probably skip (norm_residual < solution_tolerance)) { if (!quiet) std::cout << 'Newton iterations converged!' << std::endl; newton_converged = true; break; // out of Newton iterations } } // end nonlinear loop if (!newton_converged) { std::cout << 'Newton iterations of augmented system did not converge!' << std::endl; // Reduce ds_current, recompute the solution and parameter, and continue to next // arcstep, if there is one. ds_current *= 0.5; // Go back to previous solution and parameter value. *solution = *previous_u; *continuation_parameter = old_continuation_parameter; // Compute new predictor with smaller ds apply_predictor(); } else { // Set step convergence and break out arcstep_converged=true; break; // out of arclength reduction loop } } // end loop over arclength reductions // Check for convergence of the whole arcstep. If not converged at this // point, we have no choice but to quit. if (!arcstep_converged) { std::cout << 'Arcstep failed to converge after max number of reductions! Exiting...' << std::endl; libmesh_error(); } // Print converged solution control parameter and max value. std::cout << 'lambda_current=' << *continuation_parameter << std::endl; //std::cout << 'u_max=' << solution->max() << std::endl; // Reset old stream precision and flags. std::cout.precision(old_precision); std::cout.unsetf(std::ios_base::scientific); // Note: we don't want to go on to the next guess yet, since the user may // want to post-process this data. It's up to the user to call advance_arcstep() // when they are ready to go on. }
 

Number System::current_solution (const unsigned intglobal_dof_number) const [inherited]Returns:

the current solution for the specified global DOF.

Definition at line 115 of file system.C.

References System::current_local_solution, and System::n_dofs().

Referenced by ExactSolution::_compute_error(), UniformRefinementEstimator::_estimate_error(), HPCoarsenTest::add_projection(), JumpErrorEstimator::estimate_error(), ExactErrorEstimator::estimate_error(), FEMSystem::eulerian_residual(), PatchRecoveryErrorEstimator::EstimateError::operator()(), FEMContext::reinit(), HPCoarsenTest::select_refinement(), VTKIO::solution_to_vtk(), EnsightIO::write_scalar_ascii(), and EnsightIO::write_vector_ascii().

{
  // Check the sizes
  libmesh_assert (global_dof_number < _dof_map->n_dofs());
  libmesh_assert (global_dof_number < current_local_solution->size());
   
  return (*current_local_solution)(global_dof_number);
}
 

void System::deactivate () [inline, inherited]Deactivates the system. Only active systems are solved.

Definition at line 1111 of file system.h.

References System::_active.

{
  _active = false;
}
 

virtual bool DifferentiableSystem::element_constraint (boolrequest_jacobian, DiffContext &) [inline, virtual, inherited]Adds the constraint contribution on elem to elem_residual. If this method receives request_jacobian = true, then it should compute elem_jacobian and return true if possible. If elem_jacobian has not been computed then the method should return false.

Users may need to reimplement this for their particular PDE. To implement the constraint 0 = G(u), the user should examine u = elem_solution and add (G(u), phi_i) to elem_residual.

Definition at line 136 of file diff_system.h.

Referenced by SteadySolver::element_residual(), EulerSolver::element_residual(), Euler2Solver::element_residual(), and EigenTimeSolver::element_residual().

                                                  {
    return request_jacobian;
  }
 

virtual void DifferentiableSystem::element_postprocess (DiffContext &) [inline, virtual, inherited]Does any work that needs to be done on elem in a postprocessing loop.

Definition at line 204 of file diff_system.h.

Referenced by FEMSystem::postprocess().

{}
 

virtual void DifferentiableSystem::element_qoi (DiffContext &) [inline, virtual, inherited]Does any work that needs to be done on elem in a quantity of interest assembly loop, outputting to element_qoi.

Definition at line 216 of file diff_system.h.

Referenced by FEMSystem::assemble_qoi().

{}
 

virtual void DifferentiableSystem::element_qoi_derivative (DiffContext &) [inline, virtual, inherited]Does any work that needs to be done on elem in a quantity of interest derivative assembly loop, outputting to element_residual.

Definition at line 222 of file diff_system.h.

Referenced by FEMSystem::assemble_qoi_derivative().

{}
 

virtual bool DifferentiableSystem::element_time_derivative (boolrequest_jacobian, DiffContext &) [inline, virtual, inherited]Adds the time derivative contribution on elem to elem_residual. If this method receives request_jacobian = true, then it should compute elem_jacobian and return true if possible. If elem_jacobian has not been computed then the method should return false.

Users need to reimplement this for their particular PDE. To implement the physics model du/dt = F(u), the user should examine u = elem_solution and add (F(u), phi_i) to elem_residual.

Definition at line 119 of file diff_system.h.

Referenced by SteadySolver::element_residual(), EulerSolver::element_residual(), Euler2Solver::element_residual(), and EigenTimeSolver::element_residual().

                                                       {
    return request_jacobian;
  }
 

bool FEMSystem::eulerian_residual (boolrequest_jacobian, DiffContext &context) [virtual, inherited]Adds a pseudo-convection contribution on elem to elem_residual, if the nodes of elem are being translated by a moving mesh.

This function assumes that the user's time derivative equations (except for any equations involving unknown mesh xyz coordinates themselves) are expressed in an Eulerian frame of reference, and that the user is satisfied with an unstabilized convection term. Lagrangian equations will probably require overriding eulerian_residual() with a blank function; ALE or stabilized formulations will require reimplementing eulerian_residual() entirely.

Reimplemented from DifferentiableSystem.

Definition at line 846 of file fem_system.C.

References FEMSystem::_mesh_sys, FEMSystem::_mesh_x_var, FEMSystem::_mesh_y_var, FEMSystem::_mesh_z_var, DifferentiableSystem::_time_evolving, System::current_solution(), DifferentiableSystem::deltat, DiffContext::dof_indices_var, DiffContext::elem_subjacobians, DiffContext::elem_subresiduals, FEMContext::element_fe_var, FEMContext::element_qrule, AutoPtr< Tp >::get(), FEMContext::interior_gradient(), libMesh::invalid_uint, libmesh_real(), QBase::n_points(), System::n_vars(), System::number(), UnsteadySolver::old_nonlinear_solution(), and DifferentiableSystem::time_solver.

{
  // Only calculate a mesh movement residual if it's necessary
  if (_mesh_sys == libMesh::invalid_uint)
    return request_jacobian;

  FEMContext &context = libmesh_cast_ref<FEMContext&>(c);

  // This function only supports fully coupled mesh motion for now
  libmesh_assert(_mesh_sys == this->number());

  unsigned int n_qpoints = context.element_qrule->n_points();

  const unsigned int n_x_dofs = (_mesh_x_var == libMesh::invalid_uint) ?
                                0 : context.dof_indices_var[_mesh_x_var].size();
  const unsigned int n_y_dofs = (_mesh_y_var == libMesh::invalid_uint) ?
                                0 : context.dof_indices_var[_mesh_y_var].size();
  const unsigned int n_z_dofs = (_mesh_z_var == libMesh::invalid_uint) ?
                                0 : context.dof_indices_var[_mesh_z_var].size();

  const unsigned int mesh_xyz_var = n_x_dofs ? _mesh_x_var :
                                   (n_y_dofs ? _mesh_y_var :
                                   (n_z_dofs ? _mesh_z_var :
                                    libMesh::invalid_uint));

  // If we're our own _mesh_sys, we'd better be in charge of
  // at least one coordinate, and we'd better have the same
  // FE type for all coordinates we are in charge of
  libmesh_assert(mesh_xyz_var != libMesh::invalid_uint);
  libmesh_assert(!n_x_dofs || context.element_fe_var[_mesh_x_var] ==
                              context.element_fe_var[mesh_xyz_var]);
  libmesh_assert(!n_y_dofs || context.element_fe_var[_mesh_y_var] ==
                              context.element_fe_var[mesh_xyz_var]);
  libmesh_assert(!n_z_dofs || context.element_fe_var[_mesh_z_var] ==
                              context.element_fe_var[mesh_xyz_var]);

  const std::vector<std::vector<Real> >     &psi =
    context.element_fe_var[mesh_xyz_var]->get_phi();

  for (unsigned int var = 0; var != this->n_vars(); ++var)
    {
      // Mesh motion only affects time-evolving variables
      if (!_time_evolving[var])
        continue;

      // The mesh coordinate variables themselves are Lagrangian,
      // not Eulerian, and no convective term is desired.
      if (_mesh_sys == this->number() &&
          (var == _mesh_x_var ||
           var == _mesh_y_var ||
           var == _mesh_z_var))
        continue;

      // Some of this code currently relies on the assumption that
      // we can pull mesh coordinate data from our own system
      if (_mesh_sys != this->number())
        libmesh_not_implemented();

      // This residual should only be called by unsteady solvers:
      // if the mesh is steady, there's no mesh convection term!
      UnsteadySolver *unsteady =
        dynamic_cast<UnsteadySolver *>(this->time_solver.get());
      if (!unsteady)
        return request_jacobian;

      const std::vector<Real> &JxW = 
        context.element_fe_var[var]->get_JxW();

      const std::vector<std::vector<Real> >     &phi =
        context.element_fe_var[var]->get_phi();

      const std::vector<std::vector<RealGradient> > &dphi =
        context.element_fe_var[var]->get_dphi();

      const unsigned int n_u_dofs = context.dof_indices_var[var].size();

      DenseSubVector<Number> &Fu = *context.elem_subresiduals[var];
      DenseSubMatrix<Number> &Kuu = *context.elem_subjacobians[var][var];

      DenseSubMatrix<Number> *Kux = n_x_dofs ?
        context.elem_subjacobians[var][_mesh_x_var] : NULL;
      DenseSubMatrix<Number> *Kuy = n_y_dofs ?
        context.elem_subjacobians[var][_mesh_y_var] : NULL;
      DenseSubMatrix<Number> *Kuz = n_z_dofs ?
        context.elem_subjacobians[var][_mesh_z_var] : NULL;

      std::vector<Real> delta_x(n_x_dofs, 0.);
      std::vector<Real> delta_y(n_y_dofs, 0.);
      std::vector<Real> delta_z(n_z_dofs, 0.);

      for (unsigned int i = 0; i != n_x_dofs; ++i)
        {
          unsigned int j = context.dof_indices_var[_mesh_x_var][i];
          delta_x[i] = libmesh_real(this->current_solution(j)) -
                       libmesh_real(unsteady->old_nonlinear_solution(j));
        }

      for (unsigned int i = 0; i != n_y_dofs; ++i)
        {
          unsigned int j = context.dof_indices_var[_mesh_y_var][i];
          delta_y[i] = libmesh_real(this->current_solution(j)) -
                       libmesh_real(unsteady->old_nonlinear_solution(j));
        }

      for (unsigned int i = 0; i != n_z_dofs; ++i)
        {
          unsigned int j = context.dof_indices_var[_mesh_z_var][i];
          delta_z[i] = libmesh_real(this->current_solution(j)) -
                       libmesh_real(unsteady->old_nonlinear_solution(j));
        }

      for (unsigned int qp = 0; qp != n_qpoints; ++qp)
        {
          Gradient grad_u = context.interior_gradient(var, qp);
          RealGradient convection(0.);

          for (unsigned int i = 0; i != n_x_dofs; ++i)
            convection(0) += delta_x[i] * psi[i][qp];
          for (unsigned int i = 0; i != n_y_dofs; ++i)
            convection(1) += delta_y[i] * psi[i][qp];
          for (unsigned int i = 0; i != n_z_dofs; ++i)
            convection(2) += delta_z[i] * psi[i][qp];

          for (unsigned int i = 0; i != n_u_dofs; ++i)
            {
              Number JxWxPhiI = JxW[qp] * phi[i][qp];
              Fu(i) += (convection * grad_u) * JxWxPhiI;
              if (request_jacobian)
                {
                  Number JxWxPhiI = JxW[qp] * phi[i][qp];
                  for (unsigned int j = 0; j != n_u_dofs; ++j)
                    Kuu(i,j) += JxWxPhiI * (convection * dphi[j][qp]);

                  Number JxWxPhiIoverDT = JxWxPhiI/this->deltat;

                  Number JxWxPhiIxDUDXoverDT = JxWxPhiIoverDT * grad_u(0);
                  for (unsigned int j = 0; j != n_x_dofs; ++j)
                    (*Kux)(i,j) += JxWxPhiIxDUDXoverDT * psi[j][qp];

                  Number JxWxPhiIxDUDYoverDT = JxWxPhiIoverDT * grad_u(1);
                  for (unsigned int j = 0; j != n_y_dofs; ++j)
                    (*Kuy)(i,j) += JxWxPhiIxDUDYoverDT * psi[j][qp];

                  Number JxWxPhiIxDUDZoverDT = JxWxPhiIoverDT * grad_u(2);
                  for (unsigned int j = 0; j != n_z_dofs; ++j)
                    (*Kuz)(i,j) += JxWxPhiIxDUDZoverDT * psi[j][qp];
                }
            }
        }
    }

  return request_jacobian;
}
 

NumericVector< Number > & System::get_adjoint_solution () [inherited]Returns:

a reference to the system's adjoint solution vector name adjoint_solution. Adjoint system needs to be solved first.

Definition at line 637 of file system.C.

References System::get_vector().

Referenced by AdjointResidualErrorEstimator::estimate_error().

{
  // Get the adjoint solution using the get_vector function declared above
  return this->get_vector('adjoint_solution');
}
 

const NumericVector< Number > & System::get_adjoint_solution () const [inherited]Returns:

a reference to the system's adjoint solution vector name adjoint_solution. Adjoint system needs to be solved first.

Definition at line 645 of file system.C.

References System::get_vector().

{
  return this->get_vector('adjoint_solution');
}
 

const DofMap & System::get_dof_map () const [inline, inherited]Returns:

a constant reference to this system's _dof_map.

Definition at line 1079 of file system.h.

References System::_dof_map.

Referenced by __libmesh_petsc_diff_solver_jacobian(), __libmesh_petsc_diff_solver_residual(), ExactSolution::_compute_error(), UniformRefinementEstimator::_estimate_error(), HPCoarsenTest::add_projection(), PetscDiffSolver::adjoint_solve(), NewtonSolver::adjoint_solve(), FEMSystem::assemble_qoi_derivative(), FEMSystem::assembly(), EquationSystems::build_discontinuous_solution_vector(), EquationSystems::build_solution_vector(), System::calculate_norm(), DofMap::enforce_constraints_exactly(), JumpErrorEstimator::estimate_error(), ExactErrorEstimator::estimate_error(), System::get_info(), EigenSystem::init_data(), ImplicitSystem::init_matrices(), System::local_dof_indices(), DofMap::max_constraint_error(), FEMSystem::mesh_position_get(), UnsteadySolver::old_nonlinear_solution(), System::ProjectVector::operator()(), PatchRecoveryErrorEstimator::EstimateError::operator()(), ErrorVector::plot_error(), System::project_vector(), FEMContext::reinit(), EquationSystems::reinit(), HPCoarsenTest::select_refinement(), UnsteadySolver::solve(), PetscDiffSolver::solve(), NewtonSolver::solve(), EnsightIO::write_scalar_ascii(), and EnsightIO::write_vector_ascii().

{
  return *_dof_map;
}
 

DofMap & System::get_dof_map () [inline, inherited]Returns:

a writeable reference to this system's _dof_map.

Definition at line 1087 of file system.h.

References System::_dof_map.

{
  return *_dof_map;
}
 

const EquationSystems& System::get_equation_systems () const [inline, inherited]Returns:

a constant reference to this system's parent EquationSystems object.

Definition at line 252 of file system.h.

References System::_equation_systems.

Referenced by UniformRefinementEstimator::_estimate_error(), LinearImplicitSystem::adjoint_solve(), NewmarkSystem::clear(), FrequencySystem::clear_all(), ExactErrorEstimator::find_squared_element_error(), FrequencySystem::init_data(), FrequencySystem::n_frequencies(), FrequencySystem::set_current_frequency(), FrequencySystem::set_frequencies(), FrequencySystem::set_frequencies_by_range(), FrequencySystem::set_frequencies_by_steps(), NewmarkSystem::set_newmark_parameters(), NonlinearImplicitSystem::set_solver_parameters(), LinearImplicitSystem::solve(), FrequencySystem::solve(), and EigenSystem::solve().

{ return _equation_systems; }
 

EquationSystems& System::get_equation_systems () [inline, inherited]Returns:

a reference to this system's parent EquationSystems object.

Definition at line 257 of file system.h.

References System::_equation_systems.

{ return _equation_systems; }
 

std::string System::get_info () const [inherited]Returns:

a string containing information about the system.

Definition at line 1001 of file system.C.

References Utility::enum_to_string< FEFamily >(), Utility::enum_to_string< InfMapType >(), Utility::enum_to_string< Order >(), FEType::family, System::get_dof_map(), FEType::inf_map, System::n_constrained_dofs(), System::n_dofs(), System::n_local_dofs(), System::n_vars(), System::n_vectors(), System::name(), FEType::order, FEType::radial_family, FEType::radial_order, System::system_type(), System::variable_name(), and DofMap::variable_type().

{
  std::ostringstream out;

  
  const std::string& sys_name = this->name();
      
  out << '   System '' << sys_name << ''
      << '    Type ''  << this->system_type() << ''
      << '    Variables=';
  
  for (unsigned int vn=0; vn<this->n_vars(); vn++)
      out << ''' << this->variable_name(vn) << '' ';
     
  out << ';

  out << '    Finite Element Types=';
#ifndef LIBMESH_ENABLE_INFINITE_ELEMENTS
  for (unsigned int vn=0; vn<this->n_vars(); vn++)
    out << '''
        << Utility::enum_to_string<FEFamily>(this->get_dof_map().variable_type(vn).family)
        << '' ';  
#else
  for (unsigned int vn=0; vn<this->n_vars(); vn++)
    {
      out << '''
          << Utility::enum_to_string<FEFamily>(this->get_dof_map().variable_type(vn).family)
          << '', ''
          << Utility::enum_to_string<FEFamily>(this->get_dof_map().variable_type(vn).radial_family)
          << '' ';
    }

  out << ' << '    Infinite Element Mapping=';
  for (unsigned int vn=0; vn<this->n_vars(); vn++)
    out << '''
        << Utility::enum_to_string<InfMapType>(this->get_dof_map().variable_type(vn).inf_map)
        << '' ';
#endif      

  out << ';
      
  out << '    Approximation Orders=';
  for (unsigned int vn=0; vn<this->n_vars(); vn++)
    {
#ifndef LIBMESH_ENABLE_INFINITE_ELEMENTS
      out << '''
          << Utility::enum_to_string<Order>(this->get_dof_map().variable_type(vn).order)
          << '' ';
#else
      out << '''
          << Utility::enum_to_string<Order>(this->get_dof_map().variable_type(vn).order)
          << '', ''
          << Utility::enum_to_string<Order>(this->get_dof_map().variable_type(vn).radial_order)
          << '' ';
#endif
    }

  out << ';
      
  out << '    n_dofs()='             << this->n_dofs()             << ';
  out << '    n_local_dofs()='       << this->n_local_dofs()       << ';
#ifdef LIBMESH_ENABLE_AMR
  out << '    n_constrained_dofs()=' << this->n_constrained_dofs() << ';
#endif

  out << '    ' << 'n_vectors()='  << this->n_vectors()  << ';
//   out << '    ' << 'n_additional_matrices()=' << this->n_additional_matrices() << ';
  
  return out.str();
}
 

std::string ReferenceCounter::get_info () [static, inherited]Gets a string containing the reference information.

Definition at line 45 of file reference_counter.C.

References ReferenceCounter::_counts, and Quality::name().

Referenced by ReferenceCounter::print_info().

{
#if defined(LIBMESH_ENABLE_REFERENCE_COUNTING) && defined(DEBUG)

  std::ostringstream out;
  
  out << '
      << ' ---------------------------------------------------------------------------- 
      << '| Reference count information                                                |
      << ' ---------------------------------------------------------------------------- ;
  
  for (Counts::iterator it = _counts.begin();
       it != _counts.end(); ++it)
    {
      const std::string name(it->first);
      const unsigned int creations    = it->second.first;
      const unsigned int destructions = it->second.second;

      out << '| ' << name << ' reference count information:
          << '|  Creations:    ' << creations    << '
          << '|  Destructions: ' << destructions << ';
    }
  
  out << ' ---------------------------------------------------------------------------- ;

  return out.str();

#else

  return '';
  
#endif
}
 

SparseMatrix< Number > & ImplicitSystem::get_matrix (const std::string &mat_name) [inherited]Returns:

a writeable reference to this system's additional matrix named mat_name. None of these matrices is involved in the solution process. Access is only granted when the matrix is already properly initialized.

Definition at line 212 of file implicit_system.C.

References ImplicitSystem::_matrices.

{
  // Make sure the matrix exists
  matrices_iterator pos = _matrices.find (mat_name);
  
  if (pos == _matrices.end())
    {
      std::cerr << 'ERROR: matrix '
                << mat_name
                << ' does not exist in this system!'
                << std::endl;      
      libmesh_error();
    }
  
  return *(pos->second);
}
 

const SparseMatrix< Number > & ImplicitSystem::get_matrix (const std::string &mat_name) const [inherited]Returns:

a const reference to this system's additional matrix named mat_name. None of these matrices is involved in the solution process. Access is only granted when the matrix is already properly initialized.

Definition at line 193 of file implicit_system.C.

References ImplicitSystem::_matrices.

Referenced by PetscDiffSolver::adjoint_solve(), NewtonSolver::adjoint_solve(), LinearImplicitSystem::adjoint_solve(), NewmarkSystem::compute_matrix(), NewtonSolver::solve(), LinearImplicitSystem::solve(), EigenTimeSolver::solve(), and NewmarkSystem::update_rhs().

{
  // Make sure the matrix exists
  const_matrices_iterator pos = _matrices.find (mat_name);
  
  if (pos == _matrices.end())
    {
      std::cerr << 'ERROR: matrix '
                << mat_name
                << ' does not exist in this system!'
                << std::endl;      
      libmesh_error();
    }
  
  return *(pos->second);
}
 

const MeshBase & System::get_mesh () const [inline, inherited]Returns:

a constant reference to this systems's _mesh.

Definition at line 1063 of file system.h.

References System::_mesh.

Referenced by ExactSolution::_compute_error(), HPCoarsenTest::add_projection(), FEMSystem::assemble_qoi(), FEMSystem::assemble_qoi_derivative(), FEMSystem::assembly(), System::calculate_norm(), PatchRecoveryErrorEstimator::estimate_error(), JumpErrorEstimator::estimate_error(), ExactErrorEstimator::estimate_error(), AdjointResidualErrorEstimator::estimate_error(), System::init_data(), EigenSystem::init_data(), ImplicitSystem::init_matrices(), System::local_dof_indices(), DofMap::max_constraint_error(), FEMSystem::mesh_position_get(), FEMSystem::mesh_position_set(), System::ProjectVector::operator()(), PatchRecoveryErrorEstimator::EstimateError::operator()(), FEMSystem::postprocess(), System::project_vector(), System::read_header(), System::read_legacy_data(), System::read_parallel_data(), System::read_serialized_vector(), HPSingularity::select_refinement(), HPCoarsenTest::select_refinement(), System::write_header(), System::write_parallel_data(), System::write_serialized_vector(), and System::zero_variable().

{
  return _mesh;
}
 

MeshBase & System::get_mesh () [inline, inherited]Returns:

a reference to this systems's _mesh.

Definition at line 1071 of file system.h.

References System::_mesh.

{
  return _mesh;
}
 

const NumericVector< Number > & System::get_vector (const std::string &vec_name) const [inherited]Returns:

a const reference to this system's additional vector named vec_name. Access is only granted when the vector is already properly initialized.

Definition at line 551 of file system.C.

References System::_vectors.

Referenced by UniformRefinementEstimator::_estimate_error(), UnsteadySolver::advance_timestep(), AdaptiveTimeSolver::advance_timestep(), System::compare(), UnsteadySolver::du(), System::get_adjoint_solution(), NewmarkSystem::initial_conditions(), UnsteadySolver::solve(), TwostepTimeSolver::solve(), FrequencySystem::solve(), NewmarkSystem::update_rhs(), and NewmarkSystem::update_u_v_a().

{
  // Make sure the vector exists
  const_vectors_iterator pos = _vectors.find(vec_name);
  
  if (pos == _vectors.end())
    {
      std::cerr << 'ERROR: vector '
                << vec_name
                << ' does not exist in this system!'
                << std::endl;      
      libmesh_error();
    }
  
  return *(pos->second);
}
 

NumericVector< Number > & System::get_vector (const std::string &vec_name) [inherited]Returns:

a writeable reference to this system's additional vector named vec_name. Access is only granted when the vector is already properly initialized.

Definition at line 570 of file system.C.

References System::_vectors.

{
  // Make sure the vector exists
  vectors_iterator pos = _vectors.find(vec_name);
  
  if (pos == _vectors.end())
    {
      std::cerr << 'ERROR: vector '
                << vec_name
                << ' does not exist in this system!'
                << std::endl;      
      libmesh_error();
    }
  
  return *(pos->second);
}
 

const NumericVector< Number > & System::get_vector (const unsigned intvec_num) const [inherited]Returns:

a const reference to this system's additional vector number vec_num (where the vectors are counted starting with 0).

Definition at line 589 of file system.C.

References System::vectors_begin(), and System::vectors_end().

{
  const_vectors_iterator v = vectors_begin();
  const_vectors_iterator v_end = vectors_end();
  unsigned int num = 0;
  while((num<vec_num) && (v!=v_end))
    {
      num++;
      ++v;
    }
  libmesh_assert(v!=v_end);
  return *(v->second);
}
 

NumericVector< Number > & System::get_vector (const unsigned intvec_num) [inherited]Returns:

a writeable reference to this system's additional vector number vec_num (where the vectors are counted starting with 0).

Definition at line 605 of file system.C.

References System::vectors_begin(), and System::vectors_end().

{
  vectors_iterator v = vectors_begin();
  vectors_iterator v_end = vectors_end();
  unsigned int num = 0;
  while((num<vec_num) && (v!=v_end))
    {
      num++;
      ++v;
    }
  libmesh_assert(v!=v_end);
  return *(v->second);
}
 

bool System::has_variable (const std::string &var) const [inherited]Returns:

true if a variable named var exists in this System

Definition at line 703 of file system.C.

References System::_variable_numbers.

Referenced by GMVIO::copy_nodal_solution().

{
  return _variable_numbers.count(var);
}
 

bool ImplicitSystem::have_matrix (const std::string &mat_name) const [inline, inherited]Returns:

true if this System has a matrix associated with the given name, false otherwise.

Definition at line 199 of file implicit_system.h.

References ImplicitSystem::_matrices.

Referenced by ImplicitSystem::add_matrix(), PetscDiffSolver::adjoint_solve(), NewtonSolver::adjoint_solve(), LinearImplicitSystem::adjoint_solve(), EigenTimeSolver::init(), NewtonSolver::solve(), and LinearImplicitSystem::solve().

{
  return (_matrices.count(mat_name));
}
 

bool System::have_vector (const std::string &vec_name) const [inline, inherited]Returns:

true if this System has a vector associated with the given name, false otherwise.

Definition at line 1173 of file system.h.

References System::_vectors.

Referenced by System::add_vector().

{
  return (_vectors.count(vec_name));
}
 

void ReferenceCounter::increment_constructor_count (const std::string &name) [inline, protected, inherited]Increments the construction counter. Should be called in the constructor of any derived class that will be reference counted.

Definition at line 149 of file reference_counter.h.

References ReferenceCounter::_counts, Quality::name(), and Threads::spin_mtx.

Referenced by ReferenceCountedObject< Value >::ReferenceCountedObject().

{
  Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
  std::pair<unsigned int, unsigned int>& p = _counts[name];

  p.first++;
}
 

void ReferenceCounter::increment_destructor_count (const std::string &name) [inline, protected, inherited]Increments the destruction counter. Should be called in the destructor of any derived class that will be reference counted.

Definition at line 167 of file reference_counter.h.

References ReferenceCounter::_counts, Quality::name(), and Threads::spin_mtx.

Referenced by ReferenceCountedObject< Value >::~ReferenceCountedObject().

{
  Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
  std::pair<unsigned int, unsigned int>& p = _counts[name];

  p.second++;
}
 

void System::init () [inherited]Initializes degrees of freedom on the current mesh. Sets the

Definition at line 156 of file system.C.

References System::init_data(), System::n_vars(), and System::user_initialization().

{ 
  // First initialize any required data
  this->init_data();
  
  //If no variables have been added to this system
  //don't do anything
  if(!this->n_vars())
    return;
 
  // Then call the user-provided intialization function
  this->user_initialization();
}
 

void FEMSystem::init_context (DiffContext &c) [virtual, inherited]

Reimplemented from DifferentiableSystem.

Definition at line 731 of file fem_system.C.

References DifferentiableSystem::_time_evolving, FEMContext::element_fe_var, DifferentiableSystem::init_context(), and System::n_vars().

Referenced by FEMSystem::mesh_position_get().

{
  Parent::init_context(c);

  FEMContext &context = libmesh_cast_ref<FEMContext&>(c);

  // Make sure we're prepared to do mass integration
  for (unsigned int var = 0; var != this->n_vars(); ++var)
    if (_time_evolving[var])
      {
        context.element_fe_var[var]->get_JxW();
        context.element_fe_var[var]->get_phi();
      }
}
 

void ContinuationSystem::init_data () [protected, virtual]Initializes the member data fields associated with the system, so that, e.g., assemble() may be used.

Reimplemented from FEMSystem.

Definition at line 84 of file continuation_system.C.

References System::add_vector(), delta_u, du_ds, FEMSystem::init_data(), previous_du_ds, previous_u, y, y_old, and z.

{
  // Add a vector which stores the tangent 'du/ds' to the system and save its pointer.
  du_ds = &(add_vector('du_ds'));

  // Add a vector which stores the tangent 'du/ds' to the system and save its pointer.
  previous_du_ds = &(add_vector('previous_du_ds'));

  // Add a vector to keep track of the previous nonlinear solution
  // at the old value of lambda.
  previous_u = &(add_vector('previous_u'));

  // Add a vector to keep track of the temporary solution 'y' of Ay=G_{
}. y = &(add_vector('y')); // Add a vector to keep track of the 'old value' of 'y' which is the solution of Ay=G_{
}. y_old = &(add_vector('y_old')); // Add a vector to keep track of the temporary solution 'z' of Az=-G. z = &(add_vector('z')); // Add a vector to keep track of the Newton update during the constrained PDE solves. delta_u = &(add_vector('delta_u')); // Call the Parent's initialization routine. Parent::init_data(); }
 

void ImplicitSystem::init_matrices () [protected, virtual, inherited]Initializes the matrices associated with this system.

Definition at line 92 of file implicit_system.C.

References ImplicitSystem::_can_add_matrices, ImplicitSystem::_matrices, DofMap::attach_matrix(), DofMap::compute_sparsity(), System::get_dof_map(), System::get_mesh(), SparseMatrix< T >::initialized(), and ImplicitSystem::matrix.

Referenced by ImplicitSystem::init_data(), and ImplicitSystem::reinit().

{
  libmesh_assert (matrix != NULL);

  // Check for quick return in case the system matrix
  // (and by extension all the matrices) has already
  // been initialized
  if (matrix->initialized())
    return;

  // Get a reference to the DofMap
  DofMap& dof_map = this->get_dof_map();
  
  // no chance to add other matrices
  _can_add_matrices = false;
  
  // Tell the matrices about the dof map, and vice versa
  for (matrices_iterator pos = _matrices.begin();
       pos != _matrices.end(); ++pos)
    {
      libmesh_assert (!pos->second->initialized());
      dof_map.attach_matrix (*(pos->second));
    }
  
  // Compute the sparsity pattern for the current
  // mesh and DOF distribution.  This also updates
  // additional matrices, 
DofMap now knows them dof_map.compute_sparsity (this->get_mesh()); // Initialize matrices for (matrices_iterator pos = _matrices.begin(); pos != _matrices.end(); ++pos) pos->second->init (); // Set the additional matrices to 0. for (matrices_iterator pos = _matrices.begin(); pos != _matrices.end(); ++pos) pos->second->zero (); }
 

void ContinuationSystem::initialize_tangent () [private]Before starting arclength continuation, we need at least 2 prior solutions (both solution and u_previous should be filled with meaningful values) And we need to initialize the tangent vector. This only needs to be called once.

Definition at line 125 of file continuation_system.C.

References NumericVector< T >::add(), NumericVector< T >::close(), continuation_parameter, dlambda_ds, ds_current, du_ds, NumericVector< T >::l2_norm(), old_continuation_parameter, previous_u, quiet, NumericVector< T >::scale(), set_Theta(), System::solution, solve_tangent(), tangent_initialized, Theta, Theta_LOCA, update_solution(), and y.

Referenced by continuation_solve().

{
  // Be sure the tangent was not already initialized.
  libmesh_assert (!tangent_initialized);
  
  // Compute delta_s_zero, the initial arclength travelled during the
  // first step.  Here we assume that previous_u and lambda_old store
  // the previous solution and control parameter.  You may need to
  // read in an old solution (or solve the non-continuation system)
  // first and call save_current_solution() before getting here.
  
  // 1.) Compute delta_s_zero as ||u|| - ||u_old|| + ...
  // Compute norms of the current and previous solutions
//   Real norm_u          = solution->l2_norm();
//   Real norm_previous_u = previous_u->l2_norm();
  
//   if (!quiet)
//     {
//       std::cout << 'norm_u=' << norm_u << std::endl;
//       std::cout << 'norm_previous_u=' << norm_previous_u << std::endl;
//     }
  
//   if (norm_u == norm_previous_u)
//     {
//       std::cerr << 'Warning, it appears u and previous_u are the '
//              << 'same, are you sure this is correct?'
//              << 'It's possible you forgot to set one or the other...'
//              << std::endl;
//     }
  
//   Real delta_s_zero = std::sqrt(
//                              (norm_u - norm_previous_u)*(norm_u - norm_previous_u) +
//                              (*continuation_parameter-old_continuation_parameter)*
//                              (*continuation_parameter-old_continuation_parameter)
//                              );

//   // 2.) Compute delta_s_zero as ||u -u_old|| + ...
//   *delta_u = *solution;
//   delta_u->add(-1., *previous_u);
//   delta_u->close();
//   Real norm_delta_u = delta_u->l2_norm();
//   Real norm_u          = solution->l2_norm();
//   Real norm_previous_u = previous_u->l2_norm();

//   // Scale norm_delta_u by the bigger of either norm_u or norm_previous_u
//   norm_delta_u /= std::max(norm_u, norm_previous_u);
  
//   if (!quiet)
//     {
//       std::cout << 'norm_u=' << norm_u << std::endl;
//       std::cout << 'norm_previous_u=' << norm_previous_u << std::endl;
//       //std::cout << 'norm_delta_u=' << norm_delta_u << std::endl;
//       std::cout << 'norm_delta_u/max(|u|,|u_old|)=' << norm_delta_u << std::endl;
//       std::cout << '|norm_u-norm_previous_u|=' << std::abs(norm_u - norm_previous_u) << std::endl;
//     }
  
//   const Real dlambda = *continuation_parameter-old_continuation_parameter;

//   if (!quiet)
//     std::cout << 'dlambda=' << dlambda << std::endl;
  
//   Real delta_s_zero = std::sqrt(
//                              (norm_delta_u*norm_delta_u) +
//                              (dlambda*dlambda)
//                              );
  
//   if (!quiet)
//     std::cout << 'delta_s_zero=' << delta_s_zero << std::endl;

  // 1.) + 2.)
//   // Now approximate the initial tangent d(lambda)/ds
//   this->dlambda_ds = (*continuation_parameter-old_continuation_parameter) / delta_s_zero;


//   // We can also approximate the deriv. wrt s by finite differences:
//   // du/ds = (u1 - u0) / delta_s_zero.
//   // FIXME: Use delta_u from above if we decide to keep that method.
//   *du_ds = *solution;
//   du_ds->add(-1., *previous_u);
//   du_ds->scale(1./delta_s_zero);
//   du_ds->close();


  // 3.) Treating (u-previous_u)/(lambda - lambda_old) as an approximation to du/d(lambda),
  // we follow the same technique as Carnes and Shadid.
//   const Real dlambda = *continuation_parameter-old_continuation_parameter;
//   libmesh_assert (dlambda > 0.);
  
//   // Use delta_u for temporary calculation of du/d(lambda)
//   *delta_u = *solution;
//   delta_u->add(-1., *previous_u);
//   delta_u->scale(1. / dlambda);
//   delta_u->close();

//   // Determine initial normalization parameter
//   const Real solution_size = std::max(solution->l2_norm(), previous_u->l2_norm());
//   if (solution_size > 1.)
//     {
//       Theta = 1./solution_size;
      
//       if (!quiet)
//      std::cout << 'Setting Normalization Parameter Theta=' << Theta << std::endl;
//     }
  
//   // Compute d(lambda)/ds
//   // The correct sign of d(lambda)/ds should be positive, since we assume that (lambda > lambda_old)
//   // but we could always double-check that as well.
//   Real norm_delta_u = delta_u->l2_norm();
//   this->dlambda_ds = 1. / std::sqrt(1. + Theta*Theta*norm_delta_u*norm_delta_u);

//   // Finally, compute du/ds = d(lambda)/ds * du/d(lambda)
//   *du_ds = *delta_u;
//   du_ds->scale(dlambda_ds);
//   du_ds->close();


  // 4.) Use normalized arclength formula to estimate delta_s_zero
//   // Determine initial normalization parameter
//   set_Theta();

//   // Compute (normalized) delta_s_zero
//   *delta_u = *solution;
//   delta_u->add(-1., *previous_u);
//   delta_u->close();
//   Real norm_delta_u = delta_u->l2_norm();
  
//   const Real dlambda = *continuation_parameter-old_continuation_parameter;

//   if (!quiet)
//     std::cout << 'dlambda=' << dlambda << std::endl;
  
//   Real delta_s_zero = std::sqrt(
//                              (Theta_LOCA*Theta_LOCA*Theta*norm_delta_u*norm_delta_u) +
//                              (dlambda*dlambda)
//                              );
//   *du_ds = *delta_u;
//   du_ds->scale(1./delta_s_zero);
//   dlambda_ds = dlambda / delta_s_zero;
  
//   if (!quiet)
//     {
//       std::cout << 'delta_s_zero=' << delta_s_zero << std::endl;
//       std::cout << 'initial d(lambda)/ds|_0 = ' << dlambda_ds << std::endl;
//       std::cout << 'initial ||du_ds||_0 = ' << du_ds->l2_norm() << std::endl;
//     }

//   // FIXME: Also store the initial finite-differenced approximation to -du/dlambda as y.
//   // We stick to the convention of storing negative y, since that is what we typically
//   // solve for anyway.
//   *y = *delta_u;
//   y->scale(-1./dlambda);
//   y->close();

  

  // 5.) Assume dlambda/ds_0 ~ 1/sqrt(2) and determine the value of Theta_LOCA which
  // will satisfy this criterion

  // Initial change in parameter
  const Real dlambda = *continuation_parameter-old_continuation_parameter;
  libmesh_assert (dlambda != 0.0);
  
  // Ideal initial value of dlambda_ds
  dlambda_ds = 1. / std::sqrt(2.);
  if (dlambda < 0.)
    dlambda_ds *= -1.;
  
  // This also implies the initial value of ds
  ds_current = dlambda / dlambda_ds;

  if (!quiet)
    std::cout << 'Setting ds_current|_0=' << ds_current << std::endl;
  
  // Set y = -du/dlambda using finite difference approximation
  *y = *solution;
  y->add(-1., *previous_u);
  y->scale(-1./dlambda);
  y->close();
  const Real ynorm=y->l2_norm();
  
  // Finally, set the value of du_ds to be used in the upcoming
  // tangent calculation. du/ds = du/dlambda * dlambda/ds
  *du_ds = *y;
  du_ds->scale(-dlambda_ds);
  du_ds->close();

  // Determine additional solution normalization parameter
  // (Since we just set du/ds, it will be:  ||du||*||du/ds||)
  set_Theta();

  // The value of Theta_LOCA which makes dlambda_ds = 1/sqrt(2),
  // assuming our Theta = ||du||^2.
  // Theta_LOCA = std::abs(dlambda);

  // Assuming general Theta
  Theta_LOCA = std::sqrt(1./Theta/ynorm/ynorm);
  
  
  if (!quiet)
    {
      std::cout << 'Setting initial Theta_LOCA = ' << Theta_LOCA << std::endl;
      std::cout << 'Theta_LOCA^2*Theta         = ' << Theta_LOCA*Theta_LOCA*Theta << std::endl;
      std::cout << 'initial d(lambda)/ds|_0    = ' << dlambda_ds << std::endl;
      std::cout << 'initial ||du_ds||_0        = ' << du_ds->l2_norm() << std::endl;
    }


  
  // OK, we estimated the tangent at point u0.
  // Now, to estimate the tangent at point u1, we call the solve_tangent routine.

  // Set the flag which tells us the method has been initialized.
  tangent_initialized = true;

  solve_tangent();
  
  // Advance the solution and the parameter to the next value.
  update_solution();
}
 

void System::local_dof_indices (const unsigned intvar, std::set< unsigned int > &var_indices) const [inherited]Fills the std::set with the degrees of freedom on the local processor corresponding the the variable number passed in.

Definition at line 730 of file system.C.

References MeshBase::active_local_elements_begin(), MeshBase::active_local_elements_end(), DofMap::dof_indices(), DofMap::end_dof(), DofMap::first_dof(), System::get_dof_map(), and System::get_mesh().

Referenced by System::discrete_var_norm().

{
  // Make sure the set is clear
  var_indices.clear();

  std::vector<unsigned int> dof_indices;
  
  // Begin the loop over the elements
  MeshBase::const_element_iterator       el     =
    this->get_mesh().active_local_elements_begin();
  const MeshBase::const_element_iterator end_el =
    this->get_mesh().active_local_elements_end();

  const unsigned int 
    first_local = this->get_dof_map().first_dof(),
    end_local   = this->get_dof_map().end_dof();

  for ( ; el != end_el; ++el)
    {
      const Elem* elem = *el;      
      this->get_dof_map().dof_indices (elem, dof_indices, var);

      for(unsigned int i=0; i<dof_indices.size(); i++)
        {
          unsigned int dof = dof_indices[i];

          //If the dof is owned by the local processor
          if(first_local <= dof && dof < end_local)
            var_indices.insert(dof_indices[i]);
        }
    }
}
 

bool FEMSystem::mass_residual (boolrequest_jacobian, DiffContext &context) [virtual, inherited]Adds a mass vector contribution on elem to elem_residual. If this method receives request_jacobian = true, then it should compute elem_jacobian and return true if possible. If elem_jacobian has not been computed then the method should return false.

Most problems can use the FEMSystem::mass_residual implementation, which calculates the residual (u, phi_i) and jacobian (phi_i, phi_j); few users will need to reimplement this themselves. Using a custom mass matrix (e.g. for divergence-free elements or mass lumping) requires reimplementing mass_residual().

Reimplemented from DifferentiableSystem.

Definition at line 1003 of file fem_system.C.

References DifferentiableSystem::_time_evolving, DiffContext::dof_indices_var, DiffContext::elem_solution_derivative, DiffContext::elem_subjacobians, DiffContext::elem_subresiduals, FEMContext::element_fe_var, FEMContext::element_qrule, FEMContext::interior_value(), System::n_dofs(), QBase::n_points(), and System::n_vars().

{
  FEMContext &context = libmesh_cast_ref<FEMContext&>(c);

  unsigned int n_qpoints = context.element_qrule->n_points();

  for (unsigned int var = 0; var != this->n_vars(); ++var)
    {
      if (!_time_evolving[var])
        continue;

      const std::vector<Real> &JxW = 
        context.element_fe_var[var]->get_JxW();

      const std::vector<std::vector<Real> > &phi =
        context.element_fe_var[var]->get_phi();

      const unsigned int n_dofs = context.dof_indices_var[var].size();

      DenseSubVector<Number> &Fu = *context.elem_subresiduals[var];
      DenseSubMatrix<Number> &Kuu = *context.elem_subjacobians[var][var];

      for (unsigned int qp = 0; qp != n_qpoints; ++qp)
        {
          Number u = context.interior_value(var, qp);
          Number JxWxU = JxW[qp] * u;
          for (unsigned int i = 0; i != n_dofs; ++i)
            {
              Fu(i) += JxWxU * phi[i][qp];
              if (request_jacobian && context.elem_solution_derivative)
                {
                  libmesh_assert (context.elem_solution_derivative == 1.0);

                  Number JxWxPhiI = JxW[qp] * phi[i][qp];
                  Kuu(i,i) += JxWxPhiI * phi[i][qp];
                  for (unsigned int j = i+1; j < n_dofs; ++j)
                    {
                      Number Kij = JxWxPhiI * phi[j][qp];
                      Kuu(i,j) += Kij;
                      Kuu(j,i) += Kij;
                    }
                }
            }
        }
    }

  return request_jacobian;
}
 

void FEMSystem::mesh_position_get () [inherited]Tells the FEMSystem to set the degree of freedom coefficients which should correspond to mesh nodal coordinates.

Definition at line 792 of file fem_system.C.

References FEMSystem::_mesh_sys, FEMSystem::_mesh_x_var, FEMSystem::_mesh_y_var, FEMSystem::_mesh_z_var, MeshBase::active_local_elements_begin(), MeshBase::active_local_elements_end(), FEMSystem::build_context(), DofMap::dof_indices(), DiffContext::dof_indices_var, FEMContext::elem, FEMContext::elem_position_get(), DiffContext::elem_subsolutions, System::get_dof_map(), System::get_mesh(), FEMSystem::init_context(), libMesh::invalid_uint, System::n_vars(), System::number(), System::solution, and System::update().

{
  // This function makes no sense unless we've already picked out some
  // variable(s) to reflect mesh position coordinates
  if (_mesh_sys == libMesh::invalid_uint)
    libmesh_error();

  // We currently assume mesh variables are in our own system
  if (_mesh_sys != this->number())
    libmesh_not_implemented();

  // Loop over every active mesh element on this processor
  const MeshBase& mesh = this->get_mesh();

  const DofMap& dof_map = this->get_dof_map();

  MeshBase::const_element_iterator el =
    mesh.active_local_elements_begin();
  const MeshBase::const_element_iterator end_el =
    mesh.active_local_elements_end();

  AutoPtr<DiffContext> con = this->build_context();
  FEMContext &_femcontext = libmesh_cast_ref<FEMContext&>(*con);
  this->init_context(_femcontext);

  // Get the solution's mesh variables from every element
  for ( ; el != end_el; ++el)
    {
      _femcontext.elem = *el;

      // Initialize the per-variable data for elem.
      for (unsigned int i=0; i != this->n_vars(); ++i)
        {
          dof_map.dof_indices (_femcontext.elem,
                               _femcontext.dof_indices_var[i], i);
        }

      _femcontext.elem_position_get();
      if (_mesh_x_var != libMesh::invalid_uint)
        this->solution->insert(*_femcontext.elem_subsolutions[_mesh_x_var],
                               _femcontext.dof_indices_var[_mesh_x_var]);
      if (_mesh_y_var != libMesh::invalid_uint)
        this->solution->insert(*_femcontext.elem_subsolutions[_mesh_y_var],
                               _femcontext.dof_indices_var[_mesh_y_var]);
      if (_mesh_z_var != libMesh::invalid_uint)
        this->solution->insert(*_femcontext.elem_subsolutions[_mesh_z_var],
                               _femcontext.dof_indices_var[_mesh_z_var]);
    }

  // And make sure the current_local_solution is up to date too
  this->update();
}
 

void FEMSystem::mesh_position_set () [inherited]Tells the FEMSystem to set the mesh nodal coordinates which should correspond to degree of freedom coefficients.

Definition at line 392 of file fem_system.C.

References FEMSystem::_mesh_sys, MeshBase::active_elements_begin(), MeshBase::active_elements_end(), FEMSystem::build_context(), FEMContext::elem, FEMContext::elem_position_set(), System::get_mesh(), Elem::has_children(), MeshBase::is_serial(), System::number(), and FEMContext::reinit().

Referenced by FEMSystem::solve().

{
  // If we don't need to move the mesh, we're done
  if (_mesh_sys != this->number())
    return;

  const MeshBase& mesh = this->get_mesh();

  // This code won't work on a parallelized mesh yet -
  // it won't get ancestor elements right.
  libmesh_assert(mesh.is_serial());
  
  AutoPtr<DiffContext> con = this->build_context();
  FEMContext &_femcontext = libmesh_cast_ref<FEMContext&>(*con);

  // Move every mesh element we can
  MeshBase::const_element_iterator el =
    mesh.active_elements_begin();
  const MeshBase::const_element_iterator end_el =
    mesh.active_elements_end();

  for ( ; el != end_el; ++el)
    {
      _femcontext.reinit(*this, *el);

      // This code won't handle moving subactive elements
      libmesh_assert(!_femcontext.elem->has_children());

      _femcontext.elem_position_set(0.);
    }
}
 

void FEMSystem::mesh_x_position (unsigned intsysnum, unsigned intvar) [virtual, inherited]Tells the FEMSystem that variable var from system number sysnum should be used to update the x coordinate of mesh nodes, in problems where the mesh itself is expected to move in time.

The system with mesh coordinate data (which may be this system itself, for fully coupled moving mesh problems) is currently assumed to have new (end of time step) mesh coordinates stored in solution, old (beginning of time step) mesh coordinates stored in _old_nonlinear_solution, and constant velocity motion during each time step.

Activating this function ensures that local (but not neighbor!) element geometry is correctly repositioned when evaluating element residuals.

Definition at line 756 of file fem_system.C.

References FEMSystem::_mesh_sys, FEMSystem::_mesh_x_var, libMesh::invalid_uint, and System::number().

{
  if (_mesh_sys != libMesh::invalid_uint && _mesh_sys != sysnum)
    libmesh_error();
  if (sysnum != this->number())
    libmesh_not_implemented();
  _mesh_sys = sysnum;
  _mesh_x_var = var;
}
 

void FEMSystem::mesh_y_position (unsigned intsysnum, unsigned intvar) [virtual, inherited]Tells the FEMSystem that variable var from system number sysnum should be used to update the y coordinate of mesh nodes.

Definition at line 768 of file fem_system.C.

References FEMSystem::_mesh_sys, FEMSystem::_mesh_y_var, libMesh::invalid_uint, and System::number().

{
  if (_mesh_sys != libMesh::invalid_uint && _mesh_sys != sysnum)
    libmesh_error();
  if (sysnum != this->number())
    libmesh_not_implemented();
  _mesh_sys = sysnum;
  _mesh_y_var = var;
}
 

void FEMSystem::mesh_z_position (unsigned intsysnum, unsigned intvar) [virtual, inherited]Tells the FEMSystem that variable var from system number sysnum should be used to update the z coordinate of mesh nodes.

Definition at line 780 of file fem_system.C.

References FEMSystem::_mesh_sys, FEMSystem::_mesh_z_var, libMesh::invalid_uint, and System::number().

{
  if (_mesh_sys != libMesh::invalid_uint && _mesh_sys != sysnum)
    libmesh_error();
  if (sysnum != this->number())
    libmesh_not_implemented();
  _mesh_sys = sysnum;
  _mesh_z_var = var;
}
 

unsigned int System::n_active_dofs () const [inline, inherited]Returns the number of active degrees of freedom for this System.

Definition at line 1165 of file system.h.

References System::n_constrained_dofs(), and System::n_dofs().

{
  return this->n_dofs() - this->n_constrained_dofs();
}
 

unsigned int System::n_constrained_dofs () const [inherited]Returns:

the number of constrained degrees of freedom in the system.

Definition at line 93 of file system.C.

References System::_dof_map.

Referenced by System::get_info(), and System::n_active_dofs().

{
#ifdef LIBMESH_ENABLE_AMR

  return _dof_map->n_constrained_dofs();

#else

  return 0;

#endif
}
 

unsigned int System::n_dofs () const [inherited]Returns:

the number of degrees of freedom in the system

Definition at line 86 of file system.C.

References System::_dof_map.

Referenced by System::add_vector(), System::current_solution(), System::get_info(), System::init_data(), FEMSystem::mass_residual(), System::n_active_dofs(), System::ProjectVector::operator()(), System::project_vector(), System::read_legacy_data(), System::reinit(), System::restrict_vectors(), and UnsteadySolver::solve().

{
  return _dof_map->n_dofs();
}
 

unsigned int System::n_local_dofs () const [inherited]Returns:

the number of degrees of freedom local to this processor

Definition at line 108 of file system.C.

References System::_dof_map, and libMesh::processor_id().

Referenced by System::add_vector(), System::get_info(), System::init_data(), System::project_vector(), System::read_serialized_blocked_dof_objects(), System::reinit(), System::restrict_vectors(), and UnsteadySolver::solve().

{
  return _dof_map->n_dofs_on_processor (libMesh::processor_id());
}
 

unsigned int ImplicitSystem::n_matrices () const [inline, inherited]Returns:

the number of matrices handled by this system

Definition at line 206 of file implicit_system.h.

References ImplicitSystem::_matrices.

{
 return _matrices.size(); 
}
 

static unsigned int ReferenceCounter::n_objects () [inline, static, inherited]Prints the number of outstanding (created, but not yet destroyed) objects.

Definition at line 76 of file reference_counter.h.

References ReferenceCounter::_n_objects.

Referenced by System::read_serialized_blocked_dof_objects(), and System::write_serialized_blocked_dof_objects().

  { return _n_objects; }
 

unsigned int System::n_vars () const [inline, inherited]Returns:

the number of variables in the system

Definition at line 1119 of file system.h.

References System::_variables.

Referenced by UniformRefinementEstimator::_estimate_error(), System::add_variable(), EquationSystems::build_discontinuous_solution_vector(), EquationSystems::build_solution_vector(), System::calculate_norm(), DiffContext::DiffContext(), JumpErrorEstimator::estimate_error(), ExactErrorEstimator::estimate_error(), ErrorEstimator::estimate_errors(), FEMSystem::eulerian_residual(), ExactSolution::ExactSolution(), FEMContext::FEMContext(), System::get_info(), System::init(), FEMSystem::init_context(), DifferentiableSystem::init_data(), FEMSystem::mass_residual(), FEMSystem::mesh_position_get(), System::ProjectVector::operator()(), PatchRecoveryErrorEstimator::EstimateError::operator()(), System::project_vector(), System::re_update(), System::read_header(), System::read_legacy_data(), System::read_parallel_data(), System::read_serialized_vector(), System::reinit(), FEMContext::reinit(), EquationSystems::reinit(), HPCoarsenTest::select_refinement(), VTKIO::solution_to_vtk(), TecplotIO::write_ascii(), TecplotIO::write_binary(), System::write_header(), System::write_parallel_data(), System::write_serialized_vector(), and System::zero_variable().

{
  return _variables.size();
}
 

unsigned int System::n_vectors () const [inline, inherited]Returns:

the number of vectors (in addition to the solution) handled by this system This is the size of the _vectors map

Definition at line 1181 of file system.h.

References System::_vectors.

Referenced by ExplicitSystem::add_system_rhs(), System::compare(), System::get_info(), System::read_header(), and System::write_header().

{
  return _vectors.size(); 
}
 

const std::string & System::name () const [inline, inherited]Returns:

the system name.

Definition at line 1047 of file system.h.

References System::_sys_name.

Referenced by System::compare(), ExactErrorEstimator::estimate_error(), ExactSolution::ExactSolution(), ExactErrorEstimator::find_squared_element_error(), System::get_info(), System::ProjectVector::operator()(), System::project_vector(), FrequencySystem::solve(), System::user_assembly(), System::user_constrain(), System::user_initialization(), System::user_QOI(), System::user_QOI_derivative(), TecplotIO::write_binary(), System::write_header(), System::write_parallel_data(), and System::write_serialized_data().

{
  return _sys_name;
}
 

unsigned int System::number () const [inline, inherited]Returns:

the system number.

Definition at line 1055 of file system.h.

References System::_sys_number.

Referenced by System::add_variable(), FEMSystem::eulerian_residual(), FEMSystem::mesh_position_get(), FEMSystem::mesh_position_set(), FEMSystem::mesh_x_position(), FEMSystem::mesh_y_position(), FEMSystem::mesh_z_position(), FEMSystem::numerical_jacobian(), System::read_legacy_data(), System::read_parallel_data(), System::read_serialized_blocked_dof_objects(), HPCoarsenTest::select_refinement(), System::write_parallel_data(), System::write_serialized_blocked_dof_objects(), and System::zero_variable().

{
  return _sys_number;
}
 

void FEMSystem::numerical_elem_jacobian (FEMContext &context) [protected, inherited]Uses the results of multiple element_residual() calls to numerically differentiate the corresponding jacobian on an element.

Definition at line 706 of file fem_system.C.

References TimeSolver::element_residual(), and FEMSystem::numerical_jacobian().

Referenced by FEMSystem::assembly().

{
  START_LOG('numerical_elem_jacobian()', 'FEMSystem');
  this->numerical_jacobian(&TimeSolver::element_residual, context);
  STOP_LOG('numerical_elem_jacobian()', 'FEMSystem');
}
 

void FEMSystem::numerical_jacobian (TimeSolverResPtrres, FEMContext &context) [protected, inherited]Uses the results of multiple res calls to numerically differentiate the corresponding jacobian.

Definition at line 603 of file fem_system.C.

References FEMSystem::_mesh_sys, FEMSystem::_mesh_x_var, FEMSystem::_mesh_y_var, FEMSystem::_mesh_z_var, DiffContext::dof_indices, DiffContext::dof_indices_var, FEMContext::elem, DiffContext::elem_jacobian, DiffContext::elem_residual, DiffContext::elem_solution, Elem::hmin(), libMesh::invalid_uint, libmesh_real(), System::number(), FEMSystem::numerical_jacobian_h, Elem::point(), and DenseVector< T >::zero().

Referenced by FEMSystem::numerical_elem_jacobian(), and FEMSystem::numerical_side_jacobian().

{
  // Logging is done by numerical_elem_jacobian
  // or numerical_side_jacobian

  DenseVector<Number> original_residual(context.elem_residual);
  DenseVector<Number> backwards_residual(context.elem_residual);
  DenseMatrix<Number> numerical_jacobian(context.elem_jacobian);
#ifdef DEBUG
  DenseMatrix<Number> old_jacobian(context.elem_jacobian);
#endif

  Real numerical_point_h = 0.;
  if (_mesh_sys == this->number())
    numerical_point_h = numerical_jacobian_h * context.elem->hmin();

  for (unsigned int j = 0; j != context.dof_indices.size(); ++j)
    {
      // Take the 'minus' side of a central differenced first derivative
      Number original_solution = context.elem_solution(j);
      context.elem_solution(j) -= numerical_jacobian_h;

      // Make sure to catch any moving mesh terms
      // FIXME - this could be less ugly
      Real *coord = NULL;
      if (_mesh_sys == this->number())
        {
          if (_mesh_x_var != libMesh::invalid_uint)
            for (unsigned int k = 0;
                 k != context.dof_indices_var[_mesh_x_var].size(); ++k)
              if (context.dof_indices_var[_mesh_x_var][k] ==
                  context.dof_indices[j])
                coord = &(context.elem->point(k)(0));
          if (_mesh_y_var != libMesh::invalid_uint)
            for (unsigned int k = 0;
                 k != context.dof_indices_var[_mesh_y_var].size(); ++k)
              if (context.dof_indices_var[_mesh_y_var][k] == 
                  context.dof_indices[j])
                coord = &(context.elem->point(k)(1));
          if (_mesh_z_var != libMesh::invalid_uint)
            for (unsigned int k = 0;
                 k != context.dof_indices_var[_mesh_z_var].size(); ++k)
              if (context.dof_indices_var[_mesh_z_var][k] ==
                  context.dof_indices[j])
                coord = &(context.elem->point(k)(2));
        }
      if (coord)
        {
          // We have enough information to scale the perturbations
          // here appropriately
          context.elem_solution(j) = original_solution - numerical_point_h;
          *coord = libmesh_real(context.elem_solution(j));
        }

      context.elem_residual.zero();
      ((*time_solver).*(res))(false, context);
#ifdef DEBUG
      libmesh_assert(old_jacobian == context.elem_jacobian);
#endif
      backwards_residual = context.elem_residual;

      // Take the 'plus' side of a central differenced first derivative
      context.elem_solution(j) = original_solution + numerical_jacobian_h;
      if (coord)
        {
          context.elem_solution(j) = original_solution + numerical_point_h;
          *coord = libmesh_real(context.elem_solution(j));
        }
      context.elem_residual.zero();
      ((*time_solver).*(res))(false, context);
#ifdef DEBUG
      libmesh_assert(old_jacobian == context.elem_jacobian);
#endif

      context.elem_solution(j) = original_solution;
      if (coord)
        {
          *coord = libmesh_real(context.elem_solution(j));
          for (unsigned int i = 0; i != context.dof_indices.size(); ++i)
            {
              numerical_jacobian(i,j) =
                (context.elem_residual(i) - backwards_residual(i)) /
                2. / numerical_point_h;
            }
        }
      else
        {
          for (unsigned int i = 0; i != context.dof_indices.size(); ++i)
            {
              numerical_jacobian(i,j) =
                (context.elem_residual(i) - backwards_residual(i)) /
                2. / numerical_jacobian_h;
            }
        }
    }

  context.elem_residual = original_residual;
  context.elem_jacobian = numerical_jacobian;
}
 

void FEMSystem::numerical_side_jacobian (FEMContext &context) [protected, inherited]Uses the results of multiple side_residual() calls to numerically differentiate the corresponding jacobian on an element's side.

Definition at line 715 of file fem_system.C.

References FEMSystem::numerical_jacobian(), and TimeSolver::side_residual().

Referenced by FEMSystem::assembly().

{
  START_LOG('numerical_side_jacobian()', 'FEMSystem');
  this->numerical_jacobian(&TimeSolver::side_residual, context);
  STOP_LOG('numerical_side_jacobian()', 'FEMSystem');
}
 

void FEMSystem::postprocess () [virtual, inherited]Runs a postprocessing loop over all elements, and if postprocess_sides is true over all sides.

Implements DifferentiableSystem.

Definition at line 426 of file fem_system.C.

References MeshBase::active_local_elements_begin(), MeshBase::active_local_elements_end(), FEMSystem::build_context(), DifferentiableSystem::compute_internal_sides, FEMContext::elem, DifferentiableSystem::element_postprocess(), FEMSystem::fe_reinit_during_postprocess, System::get_mesh(), Elem::n_sides(), Elem::neighbor(), DifferentiableSystem::postprocess_sides, FEMContext::reinit(), FEMContext::side, FEMContext::side_fe, DifferentiableSystem::side_postprocess(), and System::update().

{
  START_LOG('postprocess()', 'FEMSystem');

  const MeshBase& mesh = this->get_mesh();

  this->update();

  AutoPtr<DiffContext> con = this->build_context();
  FEMContext &_femcontext = libmesh_cast_ref<FEMContext&>(*con);

  // Loop over every active mesh element on this processor
  MeshBase::const_element_iterator el =
    mesh.active_local_elements_begin();
  const MeshBase::const_element_iterator end_el =
    mesh.active_local_elements_end();

  for ( ; el != end_el; ++el)
    {
      _femcontext.reinit(*this, *el);

      // Optionally initialize all the interior FE objects on elem.
      // if (fe_reinit_during_postprocess)
      // _femcontext.elem_fe_reinit();
      
      this->element_postprocess(_femcontext);

      for (_femcontext.side = 0;
           _femcontext.side != _femcontext.elem->n_sides();
           ++_femcontext.side)
        {
          // Don't compute on non-boundary sides unless requested
          if (!postprocess_sides ||
              (!compute_internal_sides &&
               _femcontext.elem->neighbor(_femcontext.side) != NULL))
            continue;

          // Optionally initialize all the interior FE objects on elem/side.
          // Logging of FE::reinit is done in the FE functions
          if (fe_reinit_during_postprocess)
            {
              std::map<FEType, FEBase *>::iterator fe_end =
                _femcontext.side_fe.end();
              for (std::map<FEType, FEBase *>::iterator i =
                _femcontext.side_fe.begin();
                   i != fe_end; ++i)
                {
                  i->second->reinit(_femcontext.elem, _femcontext.side);
                }
            }

          this->side_postprocess(_femcontext);
        }
    }

  STOP_LOG('postprocess()', 'FEMSystem');
}
 

void ReferenceCounter::print_info () [static, inherited]Prints the reference information to std::cout.

Definition at line 83 of file reference_counter.C.

References ReferenceCounter::get_info().

{
#if defined(LIBMESH_ENABLE_REFERENCE_COUNTING) && defined(DEBUG)
  
  std::cout << ReferenceCounter::get_info();
  
#endif
}
 

void System::project_solution (Number fptrconst Point &p,const Parameters &parameters,const std::string &sys_name,const std::string &unknown_name, Gradient gptrconst Point &p,const Parameters &parameters,const std::string &sys_name,const std::string &unknown_name, Parameters &parameters) const [inherited]Projects the continuous functions onto the current solution.

This method projects an analytic function onto the solution via L2 projections and nodal interpolations on each element.

Definition at line 237 of file system_projection.C.

References System::current_local_solution, System::project_vector(), and System::solution.

{
  this->project_vector(fptr, gptr, parameters, *solution);

  solution->localize(*current_local_solution);
}
 

bool& System::project_solution_on_reinit (void) [inline, inherited]Tells the System whether or not to project the solution vector onto new grids when the system is reinitialized. The solution will be projected unless project_solution_on_reinit() = false is called.

Definition at line 319 of file system.h.

References System::_solution_projection.

Referenced by UniformRefinementEstimator::_estimate_error().

    { return _solution_projection; }
 

void System::project_vector (Number fptrconst Point &p,const Parameters &parameters,const std::string &sys_name,const std::string &unknown_name, Gradient gptrconst Point &p,const Parameters &parameters,const std::string &sys_name,const std::string &unknown_name, Parameters &parameters, NumericVector< Number > &new_vector) const [inherited]Projects the continuous functions onto the current mesh.

This method projects an analytic function via L2 projections and nodal interpolations on each element.

Definition at line 258 of file system_projection.C.

References NumericVector< T >::close(), DofMap::enforce_constraints_exactly(), FEType::family, System::get_dof_map(), System::get_mesh(), libMesh::n_processors(), System::n_vars(), System::name(), Threads::parallel_for(), libMesh::processor_id(), libMeshEnums::SCALAR, DofMap::SCALAR_dof_indices(), NumericVector< T >::set(), System::Variable::type(), System::variable(), and System::variable_name().

Referenced by System::project_solution(), System::project_vector(), and System::restrict_vectors().

{
  START_LOG ('project_vector()', 'System');

  Threads::parallel_for (ConstElemRange (this->get_mesh().active_local_elements_begin(),
                                         this->get_mesh().active_local_elements_end(),
                                         1000),
                         ProjectSolution(*this,
                                         fptr,
                                         gptr,
                                         parameters,
                                         new_vector)
                         );

  // Also, load values into the SCALAR dofs
  // Note: We assume that all SCALAR dofs are on the
  // processor with highest ID
  if(libMesh::processor_id() == (libMesh::n_processors()-1))
  {
    const DofMap& dof_map = this->get_dof_map();
    for (unsigned int var=0; var<this->n_vars(); var++)
      if(this->variable(var).type().family == SCALAR)
        {
          std::vector<unsigned int> SCALAR_indices;
          dof_map.SCALAR_dof_indices (SCALAR_indices, var);
          const unsigned int n_SCALAR_dofs = SCALAR_indices.size();

          for (unsigned int i=0; i<n_SCALAR_dofs; i++)
          {
            const unsigned int index = SCALAR_indices[i];

            // We pass the point (i,0,0) to the fptr to distinguish
            // the different scalars within the SCALAR variable
            Point p_i(i,0,0);
            new_vector.set( index, fptr(p_i,
                                        parameters,
                                        this->name(),
                                        this->variable_name(var))
                          );
          }
        }
  }

  new_vector.close();

#ifdef LIBMESH_ENABLE_AMR
  this->get_dof_map().enforce_constraints_exactly(*this, &new_vector);
#endif

  STOP_LOG('project_vector()', 'System');
}
 

void System::project_vector (NumericVector< Number > &vector) const [protected, inherited]Projects the vector defined on the old mesh onto the new mesh.

Definition at line 43 of file system_projection.C.

References NumericVector< T >::clone(), and System::project_vector().

{
  // Create a copy of the vector, which currently
  // contains the old data.
  AutoPtr<NumericVector<Number> >
    old_vector (vector.clone());

  // Project the old vector to the new vector
  this->project_vector (*old_vector, vector);
}
 

void System::project_vector (const NumericVector< Number > &old_v, NumericVector< Number > &new_v) const [protected, inherited]Projects the vector defined on the old mesh onto the new mesh. The original vector is unchanged and the new vector is passed through the second argument.

This method projects the vector via L2 projections or nodal interpolations on each element.

This method projects a solution from an old mesh to a current, refined mesh. The input vector old_v gives the solution on the old mesh, while the new_v gives the solution (to be computed) on the new mesh.

Definition at line 60 of file system_projection.C.

References NumericVector< T >::clear(), NumericVector< T >::close(), DofMap::enforce_constraints_exactly(), FEType::family, AutoPtr< Tp >::get(), System::get_dof_map(), System::get_mesh(), libMeshEnums::GHOSTED, NumericVector< T >::init(), NumericVector< T >::local_size(), NumericVector< T >::localize(), System::n_dofs(), System::n_local_dofs(), libMesh::n_processors(), System::n_vars(), libMeshEnums::PARALLEL, Threads::parallel_for(), Threads::parallel_reduce(), libMesh::processor_id(), libMeshEnums::SCALAR, DofMap::SCALAR_dof_indices(), System::BuildProjectionList::send_list, libMeshEnums::SERIAL, NumericVector< T >::set(), NumericVector< T >::size(), System::Variable::type(), NumericVector< T >::type(), System::BuildProjectionList::unique(), and System::variable().

{
  START_LOG ('project_vector()', 'System');

  new_v.clear();

#ifdef LIBMESH_ENABLE_AMR

  // Resize the new vector and get a serial version.
  NumericVector<Number> *new_vector_ptr = NULL;
  AutoPtr<NumericVector<Number> > new_vector_built;
  NumericVector<Number> *local_old_vector;
  AutoPtr<NumericVector<Number> > local_old_vector_built;
  const NumericVector<Number> *old_vector_ptr = NULL;

  ConstElemRange active_local_elem_range 
    (this->get_mesh().active_local_elements_begin(),
     this->get_mesh().active_local_elements_end());

  // If the old vector was uniprocessor, make the new
  // vector uniprocessor
  if (old_v.type() == SERIAL)
    {
      new_v.init (this->n_dofs(), false, SERIAL);
      new_vector_ptr = &new_v;
      old_vector_ptr = &old_v;
    }

  // Otherwise it is a parallel, distributed vector, which
  // we need to localize.
  else if (old_v.type() == PARALLEL)
    {
      // Build a send list for efficient localization
      BuildProjectionList projection_list(*this);
      Threads::parallel_reduce (active_local_elem_range, 
                                projection_list);

      // Create a sorted, unique send_list
      projection_list.unique();

      new_v.init (this->n_dofs(), this->n_local_dofs(), false, PARALLEL);
      new_vector_built = NumericVector<Number>::build();
      local_old_vector_built = NumericVector<Number>::build();
      new_vector_ptr = new_vector_built.get();
      local_old_vector = local_old_vector_built.get();
      new_vector_ptr->init(this->n_dofs(), false, SERIAL);
      local_old_vector->init(old_v.size(), false, SERIAL);
      old_v.localize(*local_old_vector, projection_list.send_list);
      local_old_vector->close();
      old_vector_ptr = local_old_vector;
    }
  else if (old_v.type() == GHOSTED)
    {
      // Build a send list for efficient localization
      BuildProjectionList projection_list(*this);
      Threads::parallel_reduce (active_local_elem_range, 
                                projection_list);

      // Create a sorted, unique send_list
      projection_list.unique();

      new_v.init (this->n_dofs(), this->n_local_dofs(),
                  this->get_dof_map().get_send_list(), false, GHOSTED);

      local_old_vector_built = NumericVector<Number>::build();
      new_vector_ptr = &new_v;
      local_old_vector = local_old_vector_built.get();
      local_old_vector->init(old_v.size(), old_v.local_size(),
                             projection_list.send_list, false, GHOSTED);
      old_v.localize(*local_old_vector, projection_list.send_list);
      local_old_vector->close();
      old_vector_ptr = local_old_vector;
    }
  else // unknown old_v.type()
    {
      std::cerr << 'ERROR: Unknown old_v.type() == ' << old_v.type() 
                << std::endl;
      libmesh_error();
    }
  
  // Note that the above will have zeroed the new_vector.
  // Just to be sure, assert that new_vector_ptr and old_vector_ptr
  // were successfully set before trying to deref them.
  libmesh_assert(new_vector_ptr);
  libmesh_assert(old_vector_ptr);

  NumericVector<Number> &new_vector = *new_vector_ptr;
  const NumericVector<Number> &old_vector = *old_vector_ptr;
    
  Threads::parallel_for (active_local_elem_range,
                         ProjectVector(*this,
                                       old_vector,
                                       new_vector)
                         );

  // Copy the SCALAR dofs from old_vector to new_vector
  // Note: We assume that all SCALAR dofs are on the
  // processor with highest ID
  if(libMesh::processor_id() == (libMesh::n_processors()-1))
  {
    const DofMap& dof_map = this->get_dof_map();
    for (unsigned int var=0; var<this->n_vars(); var++)
      if(this->variable(var).type().family == SCALAR)
        {
          // We can just map SCALAR dofs directly across
          std::vector<unsigned int> new_SCALAR_indices, old_SCALAR_indices;
          dof_map.SCALAR_dof_indices (new_SCALAR_indices, var, false);
          dof_map.SCALAR_dof_indices (old_SCALAR_indices, var, true);
          const unsigned int new_n_dofs = new_SCALAR_indices.size();

          for (unsigned int i=0; i<new_n_dofs; i++)
          {
            new_vector.set( new_SCALAR_indices[i], old_vector(old_SCALAR_indices[i]) );
          }
        }
  }

  new_vector.close();

  // If the old vector was serial, we probably need to send our values
  // to other processors
  //
  // FIXME: I'm not sure how to make a NumericVector do that without
  // creating a temporary parallel vector to use localize! - RHS
  if (old_v.type() == SERIAL)
    {
      AutoPtr<NumericVector<Number> > dist_v = NumericVector<Number>::build();
      dist_v->init(this->n_dofs(), this->n_local_dofs(), false, PARALLEL);
      dist_v->close();
    
      for (unsigned int i=0; i!=dist_v->size(); i++)
        if (new_vector(i) != 0.0)
          dist_v->set(i, new_vector(i));

      dist_v->close();

      dist_v->localize (new_v, this->get_dof_map().get_send_list());
      new_v.close();
    }
  // If the old vector was parallel, we need to update it
  // and free the localized copies
  else if (old_v.type() == PARALLEL)
    {
      // We may have to set dof values that this processor doesn't
      // own in certain special cases, like LAGRANGE FIRST or
      // HERMITE THIRD elements on second-order meshes
      for (unsigned int i=0; i!=new_v.size(); i++)
        if (new_vector(i) != 0.0)
          new_v.set(i, new_vector(i));
      new_v.close();
    }

  this->get_dof_map().enforce_constraints_exactly(*this, &new_v);

#else

  // AMR is disabled: simply copy the vector
  new_v = old_v;
  
#endif // #ifdef LIBMESH_ENABLE_AMR

  STOP_LOG('project_vector()', 'System');
}
 

void System::prolong_vectors () [virtual, inherited]Prolong vectors after the mesh has refined

Definition at line 253 of file system.C.

References System::restrict_vectors().

Referenced by EquationSystems::reinit().

{
#ifdef LIBMESH_ENABLE_AMR
  // Currently project_vector handles both restriction and prolongation
  this->restrict_vectors();
#endif
}
 

void DifferentiableSystem::qoi_parameter_sensitivity (std::vector< Number * > &parameters, std::vector< Number > &sensitivities) [virtual, inherited]Solves for the derivative of the system's quantity of interest q with respect to each parameter p in parameters. Currently uses adjoint_solve, along with finite differenced derivatives (partial q / partial p) and (partial R / partial p).

FIXME - transient sensitivities are not yet implemented. TODO - Simultaneous sensitivity calculations for multiple QoIs are not yet implemented. Analytic options for partial derivatives are not yet implemented.

Reimplemented from ExplicitSystem.

Definition at line 120 of file diff_system.C.

{
  // Get ready to fill in senstivities:
  sensitivities.clear();
  sensitivities.resize(parameters.size(), 0);

  // An introduction to the problem:
  //
  // Residual R(u(p),p) = 0
  // partial R / partial u = J = system matrix
  //
  // This implies that:
  // d/dp(R) = 0
  // (partial R / partial p) + 
  // (partial R / partial u) * (partial u / partial p) = 0

  // We first do an adjoint solve:
  // J^T * z = (partial q / partial u)

  this->adjoint_solve();

  // We use the identities:
  // dq/dp = (partial q / partial p) + (partial q / partial u) *
  //         (partial u / partial p)
  // dq/dp = (partial q / partial p) + (J^T * z) *
  //         (partial u / partial p)
  // dq/dp = (partial q / partial p) + z * J *
  //         (partial u / partial p)
 
  // Leading to our final formula:
  // dq/dp = (partial q / partial p) - z * (partial R / partial p)

  for (unsigned int i=0; i != parameters.size(); ++i)
    {
      // We currently get partial derivatives via central differencing
      Number delta_p = 1e-6;

      // (partial q / partial p) ~= (q(p+dp)-q(p-dp))/(2*dp)
      // (partial R / partial p) ~= (rhs(p+dp) - rhs(p-dp))/(2*dp)

      Number old_parameter = *parameters[i];
      // Number old_qoi = this->qoi;

      *parameters[i] = old_parameter - delta_p;
      this->assemble_qoi();
      Number qoi_minus = this->qoi;

      this->assembly(true, false);
      this->rhs->close();
      AutoPtr<NumericVector<Number> > partialR_partialp = this->rhs->clone();
      *partialR_partialp *= -1;

      *parameters[i] = old_parameter + delta_p;
      this->assemble_qoi();
      Number qoi_plus = this->qoi;
      Number partialq_partialp = (qoi_plus - qoi_minus) / (2.*delta_p);

      this->assembly(true, false);
      this->rhs->close();
      *partialR_partialp += *this->rhs;
      *partialR_partialp /= (2.*delta_p);

      // Don't leave the parameter changed
      *parameters[i] = old_parameter;

      sensitivities[i] = partialq_partialp -
                         partialR_partialp->dot(this->get_adjoint_solution());
    }

  // All parameters have been reset.
  // Don't leave the qoi or system changed - principle of least
  // surprise.
  this->assembly(true, false);
  this->rhs->close();
  this->assemble_qoi();
}
 

void System::re_update () [virtual, inherited]Re-update the local values when the mesh has changed. This method takes the data updated by update() and makes it up-to-date on the current mesh.

Definition at line 312 of file system.C.

References System::current_local_solution, Utility::iota(), System::n_vars(), and System::solution.

{
  //const std::vector<unsigned int>& send_list = _dof_map->get_send_list ();

  // If this system is empty... don't do anything!
  if(!this->n_vars())
    return;

  // Explicitly build a send_list
  std::vector<unsigned int> send_list(solution->size());
  Utility::iota (send_list.begin(), send_list.end(), 0);
  
  // Check sizes
  libmesh_assert (current_local_solution->size()       == solution->size());
  // Not true with ghosted vectors
  // libmesh_assert (current_local_solution->local_size() == solution->size());
  libmesh_assert (!send_list.empty());
  libmesh_assert (send_list.size() <= solution->size());

  // Create current_local_solution from solution.  This will
  // put a local copy of solution into current_local_solution.
  solution->localize (*current_local_solution, send_list);
}
 

void System::read_header (Xdr &io, const std::string &version, const boolread_header = true, const boolread_additional_data = true, const boolread_legacy_format = false) [inherited]Reads the basic data header for this System.

Definition at line 78 of file system_io.C.

References System::_additional_data_written, System::_can_add_vectors, System::add_variable(), System::add_vector(), System::clear(), Xdr::data(), FEType::family, System::get_mesh(), FEType::inf_map, MeshBase::mesh_dimension(), libMeshEnums::MONOMIAL, System::n_vars(), System::n_vectors(), libMesh::on_command_line(), FEType::order, libMesh::processor_id(), FEType::radial_family, FEType::radial_order, Xdr::reading(), and libMeshEnums::XYZ.

Referenced by EquationSystems::_read_impl().

{
  // This method implements the input of a
  // System object, embedded in the output of
  // an EquationSystems<T_sys>.  This warrants some 
  // documentation.  The output file essentially
  // consists of 5 sections:
  //
  // for this system
  //
  //   5.) The number of variables in the system (unsigned int)
  //
  //   for each variable in the system
  //     
  //     6.) The name of the variable (string)
  //     
  //     7.) Combined in an FEType:
  //         - The approximation order(s) of the variable 
  //           (Order Enum, cast to int/s)
  //         - The finite element family/ies of the variable 
  //           (FEFamily Enum, cast to int/s)
  // 
  //   end variable loop
  //
  //   8.) The number of additional vectors (unsigned int),      
  //
  //     for each additional vector in the system object
  // 
  //     9.) the name of the additional vector  (string)
  //
  // end system
  libmesh_assert (io.reading());
  
  // Possibly clear data structures and start from scratch.
  if (read_header)
    this->clear ();

  // Figure out if we need to read infinite element information.
  // This will be true if the version string contains ' with infinite elements'
  const bool read_ifem_info =
    (version.rfind(' with infinite elements') < version.size()) ||
    libMesh::on_command_line ('--read_ifem_systems');
  
  
  {
    // 5.) 
    // Read the number of variables in the system
    unsigned int n_vars=0;
    if (libMesh::processor_id() == 0) io.data (n_vars);
    Parallel::broadcast(n_vars);
      
    for (unsigned int var=0; var<n_vars; var++)
      {             
        // 6.)
        // Read the name of the var-th variable
        std::string var_name;     
        if (libMesh::processor_id() == 0) io.data (var_name);
        Parallel::broadcast(var_name);
                
        // 7.)
        // Read the approximation order(s) of the var-th variable 
        int order=0;      
        if (libMesh::processor_id() == 0) io.data (order);
        Parallel::broadcast(order);
        
        
        // do the same for infinite element radial_order 
        int rad_order=0;
        if (read_ifem_info)
          {
            if (libMesh::processor_id() == 0) io.data(rad_order);
            Parallel::broadcast(rad_order);
          }


        // Read the finite element type of the var-th variable 
        int fam=0;
        if (libMesh::processor_id() == 0) io.data (fam);
        Parallel::broadcast(fam);
        FEType type;
        type.order  = static_cast<Order>(order);
        type.family = static_cast<FEFamily>(fam);

        // Check for incompatibilities.  The shape function indexing was
        // changed for the monomial and xyz finite element families to 
        // simplify extension to arbitrary p.  The consequence is that 
        // old restart files will not be read correctly.  This is expected
        // to be an unlikely occurance, but catch it anyway.
        if (read_legacy_format)
          if ((type.family == MONOMIAL || type.family == XYZ) && 
              ((type.order > 2 && this->get_mesh().mesh_dimension() == 2) ||
               (type.order > 1 && this->get_mesh().mesh_dimension() == 3)))
            {
              libmesh_here();
              std::cout << '*****************************************************************
                        << '* WARNING: reading a potentially incompatible restart file!!!   *
                        << '*  contact libmesh-users@lists.sourceforge.net for more details *
                        << '*****************************************************************'
                        << std::endl;
            }
                        
        // Read additional information for infinite elements    
        int radial_fam=0;
        int i_map=0;
        if (read_ifem_info)
          {
            if (libMesh::processor_id() == 0) io.data (radial_fam);
            Parallel::broadcast(radial_fam);
            if (libMesh::processor_id() == 0) io.data (i_map);
            Parallel::broadcast(i_map);
          }
        
#ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
        
        type.radial_order  = static_cast<Order>(rad_order);
        type.radial_family = static_cast<FEFamily>(radial_fam);
        type.inf_map       = static_cast<InfMapType>(i_map);      

#endif

        if (read_header) 
          this->add_variable (var_name, type);
      }
  }

  // 8.)  
  // Read the number of additional vectors.  
  unsigned int n_vectors=0;  
  if (libMesh::processor_id() == 0) io.data (n_vectors);
  Parallel::broadcast(n_vectors);
  
  // If n_vectors > 0, this means that write_additional_data
  // was true when this file was written.  We will need to
  // make use of this fact later.
  if (n_vectors > 0)
    this->_additional_data_written = true;  
  
  for (unsigned int vec=0; vec<n_vectors; vec++)
    {
      // 9.)
      // Read the name of the vec-th additional vector
      std::string vec_name;      
      if (libMesh::processor_id() == 0) io.data (vec_name);
      Parallel::broadcast(vec_name);
      
      if (read_additional_data)
        {
          // sanity checks
          libmesh_assert(this->_can_add_vectors);
          // Some systems may have added their own vectors already
//        libmesh_assert(this->_vectors.count(vec_name) == 0);

          this->add_vector(vec_name);
        }
    }
}
 

void System::read_legacy_data (Xdr &io, const boolread_additional_data = true) [inherited]Reads additional data, namely vectors, for this System.

Definition at line 241 of file system_io.C.

References System::_additional_data_written, System::_vectors, MeshBase::active_elements_begin(), MeshBase::active_elements_end(), Xdr::data(), System::get_mesh(), DofObject::invalid_id, System::n_dofs(), System::n_vars(), MeshBase::nodes_begin(), MeshBase::nodes_end(), System::number(), libMesh::processor_id(), Xdr::reading(), System::solution, and libMesh::zero.

{
  libmesh_deprecated();
  
  // This method implements the output of the vectors
  // contained in this System object, embedded in the 
  // output of an EquationSystems<T_sys>. 
  //
  //   10.) The global solution vector, re-ordered to be node-major 
  //       (More on this later.)                                    
  //                                                                
  //      for each additional vector in the object          
  //                                                                
  //      11.) The global additional vector, re-ordered to be       
  //           node-major (More on this later.)
  libmesh_assert (io.reading());

  // read and reordering buffers
  std::vector<Number> global_vector;
  std::vector<Number> reordered_vector;

  // 10.)
  // Read and set the solution vector
  {     
    if (libMesh::processor_id() == 0) io.data (global_vector);    
    Parallel::broadcast(global_vector);
    
    // Remember that the stored vector is node-major.
    // We need to put it into whatever application-specific
    // ordering we may have using the dof_map.
    reordered_vector.resize(global_vector.size());

    //std::cout << 'global_vector.size()=' << global_vector.size() << std::endl;
    //std::cout << 'this->n_dofs()=' << this->n_dofs() << std::endl;
    
    libmesh_assert (global_vector.size() == this->n_dofs());
        
    unsigned int cnt=0;

    const unsigned int sys     = this->number();
    const unsigned int n_vars  = this->n_vars();

    for (unsigned int var=0; var<n_vars; var++)
      { 
        // First reorder the nodal DOF values
        {
          MeshBase::node_iterator
            it  = this->get_mesh().nodes_begin(),
            end = this->get_mesh().nodes_end();
        
          for (; it != end; ++it)
            for (unsigned int index=0; index<(*it)->n_comp(sys,var); index++)
              {
                libmesh_assert ((*it)->dof_number(sys, var, index) !=
                        DofObject::invalid_id);
                
                libmesh_assert (cnt < global_vector.size());
                
                reordered_vector[(*it)->dof_number(sys, var, index)] =
                  global_vector[cnt++]; 
              }
        }
        
        // Then reorder the element DOF values
        {
          MeshBase::element_iterator
            it  = this->get_mesh().active_elements_begin(),
            end = this->get_mesh().active_elements_end();

          for (; it != end; ++it)
            for (unsigned int index=0; index<(*it)->n_comp(sys,var); index++)
              {  
                libmesh_assert ((*it)->dof_number(sys, var, index) !=
                        DofObject::invalid_id);
                
                libmesh_assert (cnt < global_vector.size());
                
                reordered_vector[(*it)->dof_number(sys, var, index)] =
                  global_vector[cnt++]; 
              }
        }
      }
            
    *(this->solution) = reordered_vector;
  }

  // For each additional vector, simply go through the list.
  // ONLY attempt to do this IF additional data was actually
  // written to the file for this system (controlled by the
  // _additional_data_written flag).  
  if (this->_additional_data_written)
    {
      std::map<std::string, NumericVector<Number>* >::iterator
        pos = this->_vectors.begin();
  
      for (; pos != this->_vectors.end(); ++pos)
        {
          // 11.)          
          // Read the values of the vec-th additional vector.
          // Prior do _not_ clear, but fill with zero, since the
          // additional vectors _have_ to have the same size
          // as the solution vector
          std::fill (global_vector.begin(), global_vector.end(), libMesh::zero);

          if (libMesh::processor_id() == 0) io.data (global_vector);
          Parallel::broadcast(global_vector);

          // If read_additional_data==true, then we will keep this vector, otherwise
          // we are going to throw it away.
          if (read_additional_data)
            {
              // Remember that the stored vector is node-major.
              // We need to put it into whatever application-specific
              // ordering we may have using the dof_map.
              std::fill (reordered_vector.begin(),
                         reordered_vector.end(),
                         libMesh::zero);
        
              reordered_vector.resize(global_vector.size());    

              libmesh_assert (global_vector.size() == this->n_dofs());
        
              unsigned int cnt=0;

              const unsigned int sys     = this->number();
              const unsigned int n_vars  = this->n_vars();
        
              for (unsigned int var=0; var<n_vars; var++)
                {
                  // First reorder the nodal DOF values
                  {
                    MeshBase::node_iterator
                      it  = this->get_mesh().nodes_begin(),
                      end = this->get_mesh().nodes_end();

                    for (; it!=end; ++it)
                      for (unsigned int index=0; index<(*it)->n_comp(sys,var); index++)
                        {
                          libmesh_assert ((*it)->dof_number(sys, var, index) !=
                                  DofObject::invalid_id);
                          
                          libmesh_assert (cnt < global_vector.size());
                          
                          reordered_vector[(*it)->dof_number(sys, var, index)] =
                            global_vector[cnt++]; 
                        }
                  }

                  // Then reorder the element DOF values
                  {
                    MeshBase::element_iterator
                      it  = this->get_mesh().active_elements_begin(),
                      end = this->get_mesh().active_elements_end();

                    for (; it!=end; ++it)
                      for (unsigned int index=0; index<(*it)->n_comp(sys,var); index++)
                        {  
                          libmesh_assert ((*it)->dof_number(sys, var, index) !=
                                  DofObject::invalid_id);
                          
                          libmesh_assert (cnt < global_vector.size());
                          
                          reordered_vector[(*it)->dof_number(sys, var, index)] =
                            global_vector[cnt++]; 
                        }
                  }
                }
            
              // use the overloaded operator=(std::vector) to assign the values
              *(pos->second) = reordered_vector;
            }
        }
    } // end if (_additional_data_written)    
}
 

void System::read_parallel_data (Xdr &io, const boolread_additional_data) [inherited]Reads additional data, namely vectors, for this System. This method may safely be called on a distributed-memory mesh. This method will read an individual file for each processor in the simulation where the local solution components for that processor are stored.

This method implements the output of the vectors contained in this System object, embedded in the output of an EquationSystems<T_sys>.

9.) The global solution vector, re-ordered to be node-major (More on this later.)

for each additional vector in the object

10.) The global additional vector, re-ordered to be node-major (More on this later.)

Note that the actual IO is handled through the Xdr class (to be renamed later?) which provides a uniform interface to both the XDR (eXternal Data Representation) interface and standard ASCII output. Thus this one section of code will read XDR or ASCII files with no changes.

Definition at line 419 of file system_io.C.

References System::_vectors, Xdr::data(), System::get_mesh(), DofObject::invalid_id, Xdr::is_open(), System::n_vars(), System::number(), Xdr::reading(), and System::solution.

{
  libmesh_assert (io.reading());
  libmesh_assert (io.is_open());

  // build the ordered nodes and element maps.
  // when writing/reading parallel files we need to iterate
  // over our nodes/elements in order of increasing global id().
  // however, this is not guaranteed to be ordering we obtain
  // by using the node_iterators/element_iterators directly.
  // so build a set, sorted by id(), that provides the ordering.
  // further, for memory economy build the set but then transfer
  // its contents to vectors, which will be sorted.
  std::vector<const DofObject*> ordered_nodes, ordered_elements;
  {    
    std::set<const DofObject*, CompareDofObjectsByID>
      ordered_nodes_set (this->get_mesh().local_nodes_begin(),
                         this->get_mesh().local_nodes_end());
      
      ordered_nodes.insert(ordered_nodes.end(),
                           ordered_nodes_set.begin(),
                           ordered_nodes_set.end());
  }
  {
    std::set<const DofObject*, CompareDofObjectsByID>
      ordered_elements_set (this->get_mesh().local_elements_begin(),
                            this->get_mesh().local_elements_end());
      
      ordered_elements.insert(ordered_elements.end(),
                              ordered_elements_set.begin(),
                              ordered_elements_set.end());
  }
  
  std::vector<Number> io_buffer;
  
  // 9.)
  //
  // Actually read the solution components
  // for the ith system to disk
  io.data(io_buffer);
  
  const unsigned int sys_num = this->number();
  const unsigned int n_vars  = this->n_vars();

  unsigned int cnt=0;
  
  // Loop over each variable and each node, and read out the value.
  for (unsigned int var=0; var<n_vars; var++)
    {
      // First read the node DOF values
      for (std::vector<const DofObject*>::const_iterator
             it = ordered_nodes.begin(); it != ordered_nodes.end(); ++it)
        for (unsigned int comp=0; comp<(*it)->n_comp(sys_num, var); comp++)
          {
            libmesh_assert ((*it)->dof_number(sys_num, var, comp) !=
                    DofObject::invalid_id);
            libmesh_assert (cnt < io_buffer.size());
            this->solution->set((*it)->dof_number(sys_num, var, comp), io_buffer[cnt++]);
          }

      // Then read the element DOF values
      for (std::vector<const DofObject*>::const_iterator
             it = ordered_elements.begin(); it != ordered_elements.end(); ++it)
        for (unsigned int comp=0; comp<(*it)->n_comp(sys_num, var); comp++)
          {
            libmesh_assert ((*it)->dof_number(sys_num, var, comp) !=
                    DofObject::invalid_id);
            libmesh_assert (cnt < io_buffer.size());
            this->solution->set((*it)->dof_number(sys_num, var, comp), io_buffer[cnt++]);
          }      
    }
  
  // Only read additional vectors if wanted  
  if (read_additional_data)
    {     
      std::map<std::string, NumericVector<Number>* >::const_iterator
        pos = _vectors.begin();
  
      for(; pos != this->_vectors.end(); ++pos)
        {
          cnt=0;
          io_buffer.clear();
          
          // 10.)
          //
          // Actually read the additional vector components
          // for the ith system to disk
          io.data(io_buffer);
          
          // Loop over each variable and each node, and read out the value.
          for (unsigned int var=0; var<n_vars; var++)
            {
              // First read the node DOF values
              for (std::vector<const DofObject*>::const_iterator
                     it = ordered_nodes.begin(); it != ordered_nodes.end(); ++it)
                for (unsigned int comp=0; comp<(*it)->n_comp(sys_num, var); comp++)
                  {
                    libmesh_assert ((*it)->dof_number(sys_num, var, comp) !=
                            DofObject::invalid_id);
                    libmesh_assert (cnt < io_buffer.size());
                    this->solution->set((*it)->dof_number(sys_num, var, comp), io_buffer[cnt++]);
                  }          
              
              // Then read the element DOF values
              for (std::vector<const DofObject*>::const_iterator
                     it = ordered_elements.begin(); it != ordered_elements.end(); ++it)
                for (unsigned int comp=0; comp<(*it)->n_comp(sys_num, var); comp++)
                  {
                    libmesh_assert ((*it)->dof_number(sys_num, var, comp) !=
                            DofObject::invalid_id);
                    libmesh_assert (cnt < io_buffer.size());
                    this->solution->set((*it)->dof_number(sys_num, var, comp), io_buffer[cnt++]);
                  }           
            }
        }
    }
}
 

void System::read_serialized_data (Xdr &io, const boolread_additional_data = true) [inherited]Reads additional data, namely vectors, for this System. This method may safely be called on a distributed-memory mesh.

Definition at line 559 of file system_io.C.

References System::_vectors, Xdr::comment(), libMesh::processor_id(), System::read_serialized_vector(), and System::solution.

{
  // This method implements the input of the vectors
  // contained in this System object, embedded in the 
  // output of an EquationSystems<T_sys>. 
  //
  //   10.) The global solution vector, re-ordered to be node-major 
  //       (More on this later.)                                    
  //                                                                
  //      for each additional vector in the object          
  //                                                                
  //      11.) The global additional vector, re-ordered to be       
  //          node-major (More on this later.)
  parallel_only();
  std::string comment;

  // 10.)
  // Read the global solution vector
  {
    this->read_serialized_vector(io, *this->solution); 

    // get the comment
    if (libMesh::processor_id() == 0)
      io.comment (comment);  
  }
  
  // 11.)
  // Only read additional vectors if wanted  
  if (read_additional_data)
    {     
      std::map<std::string, NumericVector<Number>* >::const_iterator
        pos = _vectors.begin();
  
      for(; pos != this->_vectors.end(); ++pos)
        {
          this->read_serialized_vector(io, *pos->second);

          // get the comment
          if (libMesh::processor_id() == 0)
            io.comment (comment);         
            
        }
    }
}
 

void DifferentiableSystem::reinit () [virtual, inherited]Reinitializes the member data fields associated with the system, so that, e.g., assemble() may be used.

Reimplemented from ImplicitSystem.

Definition at line 67 of file diff_system.C.

References ImplicitSystem::reinit(), and DifferentiableSystem::time_solver.

{
  Parent::reinit();

  time_solver->reinit();
}
 

void System::restrict_vectors () [virtual, inherited]Restrict vectors after the mesh has coarsened

Definition at line 217 of file system.C.

References System::_dof_map, System::_solution_projection, System::_vector_projections, System::_vectors, System::current_local_solution, libMeshEnums::GHOSTED, NumericVector< T >::init(), System::n_dofs(), System::n_local_dofs(), libMeshEnums::PARALLEL, System::project_vector(), and System::solution.

Referenced by System::prolong_vectors(), and EquationSystems::reinit().

{
#ifdef LIBMESH_ENABLE_AMR
  // Restrict the _vectors on the coarsened cells
  for (vectors_iterator pos = _vectors.begin(); pos != _vectors.end(); ++pos)
    {
      NumericVector<Number>* v = pos->second;
      
      if (_vector_projections[pos->first])
        this->project_vector (*v);
      else
        v->init (this->n_dofs(), this->n_local_dofs(), false, PARALLEL);
    }

  const std::vector<unsigned int>& send_list = _dof_map->get_send_list ();

  // Restrict the solution on the coarsened cells
  if (_solution_projection)
    this->project_vector (*solution);

#ifdef LIBMESH_ENABLE_GHOSTED
  current_local_solution->init(this->n_dofs(),
                               this->n_local_dofs(), send_list, 
                               false, GHOSTED);
#else
  current_local_solution->init(this->n_dofs());
#endif

  if (_solution_projection)
    solution->localize (*current_local_solution, send_list); 

#endif // LIBMESH_ENABLE_AMR
}
 

void ContinuationSystem::save_current_solution ()Stores the current solution and continuation parameter (as 'previous_u' and 'old_continuation_paramter') for later referral. You may need to call this e.g. after the first regular solve, in order to store the first solution, before computing a second solution and beginning arclength continuation.

Definition at line 1373 of file continuation_system.C.

References continuation_parameter, old_continuation_parameter, previous_u, and System::solution.

Referenced by update_solution().

{
  // Save the old solution vector
  *previous_u = *solution;

  // Save the old value of lambda
  old_continuation_parameter = *continuation_parameter;
}
 

void ContinuationSystem::set_max_arclength_stepsize (Realmaxds) [inline]Sets (initializes) the max-allowable ds value and the current ds value. Call this before beginning arclength continuation. The default max stepsize is 0.1

Definition at line 128 of file continuation_system.h.

References ds, and ds_current.

{ ds=maxds; ds_current=maxds; }
 

void ContinuationSystem::set_Theta () [private]A centralized function for setting the normalization parameter theta

Definition at line 1074 of file continuation_system.C.

References quiet, and Theta.

Referenced by initialize_tangent(), and update_solution().

{
  // // Use the norm of the latest solution, squared.
  //const Real normu = solution->l2_norm();
  //libmesh_assert (normu != 0.0);
  //Theta = 1./normu/normu;
  
  // // 1.) Use the norm of du, squared
//   *delta_u = *solution;
//   delta_u->add(-1, *previous_u);
//   delta_u->close();
//   const Real normdu = delta_u->l2_norm();
  
//   if (normdu < 1.) // don't divide by zero or make a huge scaling parameter.
//     Theta = 1.;
//   else
//     Theta = 1./normdu/normdu;

  // 2.) Use 1.0, i.e. don't scale
  Theta=1.;

  // 3.) Use a formula which attempts to make the 'solution triangle' isosceles.
//   libmesh_assert (std::abs(dlambda_ds) < 1.);

//   *delta_u = *solution;
//   delta_u->add(-1, *previous_u);
//   delta_u->close();
//   const Real normdu = delta_u->l2_norm();

//   Theta = std::sqrt(1. - dlambda_ds*dlambda_ds) / normdu * tau * ds;


//   // 4.) Use the norm of du and the norm of du/ds
//   *delta_u = *solution;
//   delta_u->add(-1, *previous_u);
//   delta_u->close();
//   const Real normdu   = delta_u->l2_norm();
//   du_ds->close();
//   const Real normduds = du_ds->l2_norm();

//   if (normduds < 1.e-12)
//     {
//       std::cout << 'Setting initial Theta= 1./normdu/normdu' << std::endl;
//       std::cout << 'normdu=' << normdu << std::endl;

//       // Don't use this scaling if the solution delta is already O(1)
//       if (normdu > 1.)
//      Theta = 1./normdu/normdu;
//       else
//      Theta = 1.;
//     }
//   else
//     {
//       std::cout << 'Setting Theta= 1./normdu/normduds' << std::endl;
//       std::cout << 'normdu=' << normdu << std::endl;
//       std::cout << 'normduds=' << normduds << std::endl;

//       // Don't use this scaling if the solution delta is already O(1)
//       if ((normdu>1.) || (normduds>1.))
//      Theta = 1./normdu/normduds;
//       else
//      Theta = 1.;
//     }
  
 if (!quiet)
   std::cout << 'Setting Normalization Parameter Theta=' << Theta << std::endl;
}
 

void ContinuationSystem::set_Theta_LOCA () [private]A centralized function for setting the other normalization parameter, i.e. the one suggested by the LOCA developers.

Definition at line 1144 of file continuation_system.C.

References dlambda_ds, quiet, and Theta_LOCA.

Referenced by update_solution().

{
  // We also recompute the LOCA normalization parameter based on the
  // most recently computed value of dlambda_ds
  // if (!quiet)
  //   std::cout << '(Theta_LOCA) dlambda_ds=' << dlambda_ds << std::endl;

  // Formula makes no sense if |dlambda_ds| > 1
  libmesh_assert (std::abs(dlambda_ds) < 1.);
  
  // 1.) Attempt to implement the method in LOCA paper
//   const Real g = 1./std::sqrt(2.); // 'desired' dlambda_ds

//   // According to the LOCA people, we only renormalize for
//   // when |dlambda_ds| exceeds some pre-selected maximum (which they take to be zero, btw).
//   if (std::abs(dlambda_ds) > .9)
//     {
//       // Note the *= ... This is updating the previous value of Theta_LOCA
//       // Note: The LOCA people actually use Theta_LOCA^2 to normalize their arclength constraint.
//       Theta_LOCA *= std::abs( (dlambda_ds/g)*std::sqrt( (1.-g*g) / (1.-dlambda_ds*dlambda_ds) ) );
  
//       // Suggested max-allowable value for Theta_LOCA
//       if (Theta_LOCA > 1.e8)
//      {
//        Theta_LOCA = 1.e8;

//        if (!quiet)
//          std::cout << 'max Theta_LOCA=' << Theta_LOCA << ' has been selected.' << std::endl;
//      }
//     }
//   else
//     Theta_LOCA=1.0;

  // 2.) FIXME: Should we do *= or just =?  This function is of dlambda_ds is
  //  < 1,  |dlambda_ds| < 1/sqrt(2) ~~ .7071
  //  > 1,  |dlambda_ds| > 1/sqrt(2) ~~ .7071
  Theta_LOCA *= std::abs( dlambda_ds / std::sqrt( (1.-dlambda_ds*dlambda_ds) ) );

  // Suggested max-allowable value for Theta_LOCA.  I've never come close
  // to this value in my code.
  if (Theta_LOCA > 1.e8)
    {
      Theta_LOCA = 1.e8;
      
      if (!quiet)
        std::cout << 'max Theta_LOCA=' << Theta_LOCA << ' has been selected.' << std::endl;
    }
    
  // 3.) Use 1.0, i.e. don't scale
  //Theta_LOCA=1.0;

  if (!quiet)
    std::cout << 'Setting Theta_LOCA=' << Theta_LOCA << std::endl;
}
 

virtual bool DifferentiableSystem::side_constraint (boolrequest_jacobian, DiffContext &) [inline, virtual, inherited]Adds the time derivative contribution on side of elem to elem_residual. If this method receives request_jacobian = true, then it should compute elem_jacobian and return true if possible. If elem_jacobian has not been computed then the method should return false.

Users may need to reimplement this for their particular PDE depending on the boundary conditions.

Definition at line 176 of file diff_system.h.

Referenced by SteadySolver::side_residual(), EulerSolver::side_residual(), Euler2Solver::side_residual(), and EigenTimeSolver::side_residual().

                                               {
    return request_jacobian;
  }
 

virtual bool DifferentiableSystem::side_mass_residual (boolrequest_jacobian, DiffContext &) [inline, virtual, inherited]Adds a mass vector contribution on side of elem to elem_residual. If this method receives request_jacobian = true, then it should compute elem_jacobian and return true if possible. If elem_jacobian has not been computed then the method should return false.

For most problems, the default implementation of 'do nothing' is correct; users with boundary conditions including time derivatives may need to reimplement this themselves.

Definition at line 294 of file diff_system.h.

Referenced by EulerSolver::side_residual(), Euler2Solver::side_residual(), and EigenTimeSolver::side_residual().

                                                  {
    return request_jacobian;
  }
 

virtual void DifferentiableSystem::side_postprocess (DiffContext &) [inline, virtual, inherited]Does any work that needs to be done on side of elem in a postprocessing loop.

Definition at line 210 of file diff_system.h.

Referenced by FEMSystem::postprocess().

{}
 

virtual void DifferentiableSystem::side_qoi (DiffContext &) [inline, virtual, inherited]Does any work that needs to be done on side of elem in a quantity of interest assembly loop, outputting to element_qoi.

Definition at line 228 of file diff_system.h.

Referenced by FEMSystem::assemble_qoi().

{}
 

virtual void DifferentiableSystem::side_qoi_derivative (DiffContext &) [inline, virtual, inherited]Does any work that needs to be done on side of elem in a quantity of interest derivative assembly loop, outputting to element_residual.

Definition at line 235 of file diff_system.h.

Referenced by FEMSystem::assemble_qoi_derivative().

{}
 

virtual bool DifferentiableSystem::side_time_derivative (boolrequest_jacobian, DiffContext &) [inline, virtual, inherited]Adds the time derivative contribution on side of elem to elem_residual. If this method receives request_jacobian = true, then it should compute elem_jacobian and return true if possible. If elem_jacobian has not been computed then the method should return false.

Users may need to reimplement this for their particular PDE depending on the boundary conditions.

Definition at line 160 of file diff_system.h.

Referenced by SteadySolver::side_residual(), EulerSolver::side_residual(), Euler2Solver::side_residual(), and EigenTimeSolver::side_residual().

                                                    {
    return request_jacobian;
  }
 

void ContinuationSystem::solve () [virtual]Perform a standard 'solve' of the system, without doing continuation.

Reimplemented from FEMSystem.

Definition at line 115 of file continuation_system.C.

References Residual, and rhs_mode.

{
  // Set the Residual RHS mode, and call the normal solve routine.
  rhs_mode      = Residual;
  DifferentiableSystem::solve();
}
 

void ContinuationSystem::solve_tangent () [private]Special solve algorithm for solving the tangent system.

Definition at line 965 of file continuation_system.C.

References NumericVector< T >::add(), FEMSystem::assembly(), NumericVector< T >::close(), continuation_parameter, delta_u, dlambda_ds, NumericVector< T >::dot(), du_ds, G_Lambda, AutoPtr< Tp >::get(), NumericVector< T >::l2_norm(), libmesh_real(), linear_solver, ImplicitSystem::matrix, DiffSolver::max_linear_iterations, newton_solver, old_continuation_parameter, previous_dlambda_ds, previous_du_ds, previous_u, quiet, ExplicitSystem::rhs, rhs_mode, NumericVector< T >::scale(), System::solution, tangent_initialized, Theta, Theta_LOCA, DifferentiableSystem::time_solver, y, and NumericVector< T >::zero().

Referenced by advance_arcstep(), and initialize_tangent().

{
  // We shouldn't call this unless the current tangent already makes sense.
  libmesh_assert (tangent_initialized);
  
  // Set pointer to underlying Newton solver
  if (!newton_solver)
    newton_solver =
      libmesh_cast_ptr<NewtonSolver*> (this->time_solver->diff_solver().get());

  // Assemble the system matrix AND rhs, with rhs = G_{
} this->rhs_mode = G_Lambda; // Assemble Residual and Jacobian this->assembly(true, // Residual true); // Jacobian // Not sure if this is really necessary rhs->close(); // Solve G_u*y = G_{
} std::pair<unsigned int, Real> rval = linear_solver->solve(*matrix, *y, *rhs, 1.e-12, // relative linear tolerance 2*newton_solver->max_linear_iterations); // max linear iterations // FIXME: If this doesn't converge at all, the new tangent vector is // going to be really bad... if (!quiet) std::cout << 'G_u*y = G_{lambda} solver converged at step ' << rval.first << ' linear tolerance = ' << rval.second << '.' << std::endl; // Save old solution and parameter tangents for possible use in higher-order // predictor schemes. previous_dlambda_ds = dlambda_ds; *previous_du_ds = *du_ds; // 1.) Previous, probably wrong, technique! // // Solve for the updated d(lambda)/ds // // denom = N_{lambda} - (du_ds)^t y // // = d(lambda)/ds - (du_ds)^t y // Real denom = dlambda_ds - du_ds->dot(*y); // //std::cout << 'denom=' << denom << std::endl; // libmesh_assert (denom != 0.0); // dlambda_ds = 1.0 / denom; // if (!quiet) // std::cout << 'dlambda_ds=' << dlambda_ds << std::endl; // // Compute the updated value of du/ds = -_dlambda_ds * y // du_ds->zero(); // du_ds->add(-dlambda_ds, *y); // du_ds->close(); // 2.) From Brian Carnes' paper... // According to Carnes, y comes from solving G_u * y = -G_{
} y->scale(-1.); const Real ynorm = y->l2_norm(); dlambda_ds = 1. / std::sqrt(1. + Theta_LOCA*Theta_LOCA*Theta*ynorm*ynorm); // Determine the correct sign for dlambda_ds. // We will use delta_u to temporarily compute this sign. *delta_u = *solution; delta_u->add(-1., *previous_u); delta_u->close(); const Real sgn_dlambda_ds = libmesh_real(Theta_LOCA*Theta_LOCA*Theta*y->dot(*delta_u) + (*continuation_parameter-old_continuation_parameter)); if (sgn_dlambda_ds < 0.) { if (!quiet) std::cout << 'dlambda_ds is negative.' << std::endl; dlambda_ds *= -1.; } // Finally, set the new tangent vector, du/ds = dlambda/ds * y. du_ds->zero(); du_ds->add(dlambda_ds, *y); du_ds->close(); if (!quiet) { std::cout << 'd(lambda)/ds = ' << dlambda_ds << std::endl; std::cout << '||du_ds|| = ' << du_ds->l2_norm() << std::endl; } // Our next solve expects y ~ -du/dlambda, so scale it back by -1 again now. y->scale(-1.); y->close(); }
 

sys_type& ImplicitSystem::system () [inline, inherited]Returns:

a clever pointer to the system.

Reimplemented from ExplicitSystem.

Reimplemented in LinearImplicitSystem, and NonlinearImplicitSystem.

Definition at line 71 of file implicit_system.h.

{ return *this; }
 

virtual std::string ImplicitSystem::system_type () const [inline, virtual, inherited]Assembles & solves the linear system Ax=b.

Returns:

'Implicit'. Helps in identifying the system type in an equation system file.

Reimplemented from ExplicitSystem.

Reimplemented in FrequencySystem, LinearImplicitSystem, NewmarkSystem, and NonlinearImplicitSystem.

Definition at line 106 of file implicit_system.h.

{ return 'Implicit'; }
 

void FEMSystem::time_evolving (unsigned intvar) [virtual, inherited]Tells the FEMSystem that variable var is evolving with respect to time. In general, the user's init() function should call time_evolving() for any variables which behave like du/dt = F(u), and should not call time_evolving() for any variables which behave like 0 = G(u).

Most derived systems will not have to reimplment this function; however any system which reimplements mass_residual() may have to reimplement time_evolving() to prepare data structures.

Reimplemented from DifferentiableSystem.

Definition at line 748 of file fem_system.C.

References DifferentiableSystem::time_evolving().

{
  // Call the parent function
  Parent::time_evolving(var);
}
 

void System::update () [virtual, inherited]Update the local values to reflect the solution on neighboring processors.

Definition at line 293 of file system.C.

References System::_dof_map, System::current_local_solution, and System::solution.

Referenced by __libmesh_petsc_diff_solver_jacobian(), __libmesh_petsc_diff_solver_residual(), UniformRefinementEstimator::_estimate_error(), PetscDiffSolver::adjoint_solve(), NonlinearImplicitSystem::adjoint_solve(), NewtonSolver::adjoint_solve(), LinearImplicitSystem::adjoint_solve(), FEMSystem::assemble_qoi(), FEMSystem::assemble_qoi_derivative(), FEMSystem::assembly(), Problem_Interface::computeF(), GMVIO::copy_nodal_solution(), FEMSystem::mesh_position_get(), FEMSystem::postprocess(), NonlinearImplicitSystem::solve(), NewtonSolver::solve(), LinearImplicitSystem::solve(), and ExplicitSystem::solve().

{
  const std::vector<unsigned int>& send_list = _dof_map->get_send_list ();

  // Check sizes
  libmesh_assert (current_local_solution->size() == solution->size());
// More processors than elements => empty send_list
//  libmesh_assert (!send_list.empty());
  libmesh_assert (send_list.size() <= solution->size());

  // Create current_local_solution from solution.  This will
  // put a local copy of solution into current_local_solution.
  // Only the necessary values (specified by the send_list)
  // are copied to minimize communication
  solution->localize (*current_local_solution, send_list); 
}
 

void System::update_global_solution (std::vector< Number > &global_soln) const [inherited]Fill the input vector global_soln so that it contains the global solution on all processors. Requires communication with all other processors.

Definition at line 511 of file system.C.

References System::solution.

Referenced by ExactSolution::_compute_error(), EquationSystems::build_discontinuous_solution_vector(), EquationSystems::build_solution_vector(), and ExactErrorEstimator::estimate_error().

{
  global_soln.resize (solution->size());

  solution->localize (global_soln);
}
 

void System::update_global_solution (std::vector< Number > &global_soln, const unsigned intdest_proc) const [inherited]Fill the input vector global_soln so that it contains the global solution on processor dest_proc. Requires communication with all other processors.

Definition at line 520 of file system.C.

References System::solution.

{
  global_soln.resize        (solution->size());

  solution->localize_to_one (global_soln, dest_proc);
}
 

void ContinuationSystem::update_solution () [private]This function (which assumes the most recent tangent vector has been computed) updates the solution and the control parameter with the initial guess for the next point on the continuation path.

Definition at line 1201 of file continuation_system.C.

References NumericVector< T >::add(), apply_predictor(), NumericVector< T >::close(), continuation_parameter, delta_u, NumericVector< T >::dot(), ds, ds_current, ds_min, NumericVector< T >::l2_norm(), DiffSolver::max_nonlinear_iterations, newton_solver, newton_step, newton_stepgrowth_aggressiveness, old_continuation_parameter, previous_ds, previous_u, quiet, save_current_solution(), set_Theta(), set_Theta_LOCA(), System::solution, tangent_initialized, Theta, Theta_LOCA, y, and y_old.

Referenced by advance_arcstep(), and initialize_tangent().

{
  // Set some stream formatting flags
  unsigned int old_precision = std::cout.precision();
  std::cout.precision(16);
  std::cout.setf(std::ios_base::scientific);
  
  // We must have a tangent that makes sense before we can update the solution.
  libmesh_assert (tangent_initialized);

  // Compute tau, the stepsize scaling parameter which attempts to
  // reduce ds when the angle between the most recent two tangent
  // vectors becomes large.  tau is actually the (absolute value of
  // the) cosine of the angle between these two vectors... so if tau ~
  // 0 the angle is ~ 90 degrees, while if tau ~ 1 the angle is ~ 0
  // degrees.
  y_old->close();
  y->close();
  const Real yoldnorm = y_old->l2_norm();
  const Real ynorm = y->l2_norm();
  const Number yoldy = y_old->dot(*y);
  const Real yold_over_y = yoldnorm/ynorm;
  
  if (!quiet)
    {
      std::cout << 'yoldnorm=' << yoldnorm << std::endl;
      std::cout << 'ynorm='    << ynorm << std::endl;
      std::cout << 'yoldy='    << yoldy << std::endl;
      std::cout << 'yoldnorm/ynorm=' << yoldnorm/ynorm << std::endl;
    }      

  // Save the current value of ds before updating it
  previous_ds = ds_current;
  
  // // 1.) Cosine method (for some reason this always predicts the angle is ~0)
  // // Don't try divinding by zero
  // if ((yoldnorm > 1.e-12) && (ynorm > 1.e-12))
  //   tau = std::abs(yoldy) / yoldnorm  / ynorm;
  // else
  //   tau = 1.;

  // // 2.) Relative size of old and new du/dlambda method with cutoff of 0.9
  // if ((yold_over_y < 0.9) && (yold_over_y > 1.e-6))
  //   tau = yold_over_y;
  // else
  //   tau = 1.;

  // 3.) Grow (or shrink) the arclength stepsize by the ratio of du/dlambda, but do not
  // exceed the user-specified value of ds.
  if (yold_over_y > 1.e-6)
    {
      // // 1.) Scale current ds by the ratio of successive tangents.
      //       ds_current *= yold_over_y;
      //       if (ds_current > ds)
      //        ds_current = ds;

      // 2.) Technique 1 tends to shrink the step fairly well (and even if it doesn't
      // get very small, we still have step reduction) but it seems to grow the step
      // very slowly.  Another possible technique is step-doubling:
//       if (yold_over_y > 1.)
//              ds_current *= 2.;
//       else
//      ds_current *= yold_over_y;

      // 3.) Technique 2 may over-zealous when we are also using the Newton stepgrowth
      // factor.  For technique 3 we multiply by yold_over_y unless yold_over_y > 2
      // in which case we use 2.
      //       if (yold_over_y > 2.)
      //        ds_current *= 2.;
      //       else
      //        ds_current *= yold_over_y;

      // 4.) Double-or-halve.  We double the arc-step if the ratio of successive tangents
      // is larger than 'double_threshold', halve it if it is less than 'halve_threshold'
      const Real double_threshold = 0.5;
      const Real halve_threshold  = 0.5;
      if (yold_over_y > double_threshold)
        ds_current *= 2.;
      else if (yold_over_y < halve_threshold)
        ds_current *= 0.5;
        
      
      // Also possibly use the number of Newton iterations required to compute the previous
      // step (relative to the maximum-allowed number of Newton iterations) to grow the step.
      if (newton_stepgrowth_aggressiveness > 0.)
        {
          libmesh_assert (newton_solver != NULL);
          const unsigned int Nmax = newton_solver->max_nonlinear_iterations;

          // // The LOCA Newton step growth technique (note: only grows step length)
          // const Real stepratio = static_cast<Real>(Nmax-(newton_step+1))/static_cast<Real>(Nmax-1.);
          // const Real newtonstep_growthfactor = 1. + newton_stepgrowth_aggressiveness*stepratio*stepratio;

          // The 'Nopt/N' method, may grow or shrink the step.  Assume Nopt=Nmax/2.
          const Real newtonstep_growthfactor =
            newton_stepgrowth_aggressiveness * 0.5 *
            static_cast<Real>(Nmax) / static_cast<Real>(newton_step+1);

          if (!quiet)
            std::cout << 'newtonstep_growthfactor=' << newtonstep_growthfactor << std::endl;
        
          ds_current *= newtonstep_growthfactor;
        }
    }

  
  // Don't let the stepsize get above the user's maximum-allowed stepsize.
  if (ds_current > ds)
    ds_current = ds;
  
  // Check also for a minimum allowed stepsize.
  if (ds_current < ds_min)
    {
      std::cout << 'Enforcing minimum-allowed arclength stepsize of ' << ds_min << std::endl;
      ds_current = ds_min;
    }
  
  if (!quiet)
    {
      std::cout << 'Current step size: ds_current=' << ds_current << std::endl;
    }
  
  // Recompute scaling factor Theta for
  // the current solution before updating.
  set_Theta();

  // Also, recompute the LOCA scaling factor, which attempts to
  // maintain a reasonable value of dlambda/ds
  set_Theta_LOCA();
    
  std::cout << 'Theta*Theta_LOCA^2=' << Theta*Theta_LOCA*Theta_LOCA << std::endl;

  // Based on the asymptotic singular behavior of du/dlambda near simple turning points,
  // we can compute a single parameter which may suggest that we are close to a singularity.
  *delta_u = *solution;
  delta_u->add(-1, *previous_u);
  delta_u->close();
  const Real normdu   = delta_u->l2_norm();
  const Real C = (std::log (Theta_LOCA*normdu) /
                  std::log (std::abs(*continuation_parameter-old_continuation_parameter))) - 1.0;
  if (!quiet)
    std::cout << 'C=' << C << std::endl;
  
  // Save the current value of u and lambda before updating.
  save_current_solution();

  if (!quiet)
    {
      std::cout << 'Updating the solution with the tangent guess.' << std::endl;
      std::cout << '||u_old||=' << this->solution->l2_norm() << std::endl;
      std::cout << 'lambda_old=' << *continuation_parameter << std::endl;
    }

  // Since we solved for the tangent vector, now we can compute an
  // initial guess for the new solution, and an initial guess for the
  // new value of lambda.
  apply_predictor();

  if (!quiet)
    {
      std::cout << '||u_new||=' << this->solution->l2_norm() << std::endl;
      std::cout << 'lambda_new=' << *continuation_parameter << std::endl;
    }
  
  // Unset previous stream flags
  std::cout.precision(old_precision);
  std::cout.unsetf(std::ios_base::scientific);
}
 

void System::user_assembly () [virtual, inherited]Calls user's attached assembly function, or is overloaded by the user in derived classes.

Definition at line 1133 of file system.C.

References System::_assemble_system, System::_equation_systems, and System::name().

Referenced by System::assemble().

{
  // Call the user-provided assembly function,
  // if it was provided
  if (_assemble_system != NULL)
    this->_assemble_system (_equation_systems, this->name());
}
 

void System::user_constrain () [virtual, inherited]Calls user's attached constraint function, or is overloaded by the user in derived classes.

Definition at line 1143 of file system.C.

References System::_constrain_system, System::_equation_systems, and System::name().

Referenced by System::init_data(), and EquationSystems::reinit().

{
  // Call the user-provided constraint function, 
  // if it was provided
  if(_constrain_system!= NULL)
    this->_constrain_system(_equation_systems, this->name());
}
 

void System::user_initialization () [virtual, inherited]Calls user's attached initialization function, or is overloaded by the user in derived classes.

Definition at line 1124 of file system.C.

References System::_equation_systems, System::_init_system, and System::name().

Referenced by System::init(), and NewmarkSystem::initial_conditions().

{
  // Call the user-provided intialization function,
  // if it was provided
  if (_init_system != NULL)
    this->_init_system (_equation_systems, this->name());
}
 

void System::user_QOI () [virtual, inherited]Calls user's attached quantity of interest function, or is overloaded by the user in derived classes.

Definition at line 1153 of file system.C.

References System::_equation_systems, System::_qoi_evaluate, and System::name().

Referenced by System::assemble_qoi().

{
  // Call the user-provided quantity of interest function, 
  // if it was provided
  if(_qoi_evaluate != NULL)
    this->_qoi_evaluate(_equation_systems, this->name());
}
 

void System::user_QOI_derivative () [virtual, inherited]Calls user's attached quantity of interest derivative function, or is overloaded by the user in derived classes.

Definition at line 1163 of file system.C.

References System::_equation_systems, System::_qoi_evaluate_derivative, and System::name().

Referenced by System::assemble_qoi_derivative().

{
  // Call the user-provided quantity of interest derivative, 
  // if it was provided
  if(_qoi_evaluate_derivative != NULL)
    this->_qoi_evaluate_derivative(_equation_systems, this->name());
}
 

const System::Variable & System::variable (unsigned intvar) const [inline, inherited]Return a constant reference to Variable var.

Definition at line 1127 of file system.h.

References System::_variables.

Referenced by EquationSystems::build_solution_vector(), and System::project_vector().

{
  libmesh_assert (i < _variables.size());

  return _variables[i];
}
 

const std::string & System::variable_name (const unsigned inti) const [inline, inherited]Returns:

the name of variable i.

Definition at line 1137 of file system.h.

References System::_variables.

Referenced by System::add_variable(), KellyErrorEstimator::boundary_side_integration(), DiscontinuityMeasure::boundary_side_integration(), ExactErrorEstimator::estimate_error(), ExactSolution::ExactSolution(), System::get_info(), System::ProjectVector::operator()(), System::project_vector(), VTKIO::solution_to_vtk(), and System::write_header().

{
  libmesh_assert (i < _variables.size());

  return _variables[i].name();
}
 

unsigned short int System::variable_number (const std::string &var) const [inherited]Returns:

the variable number assoicated with the user-specified variable named var.

Definition at line 710 of file system.C.

References System::_variable_numbers, and System::_variables.

Referenced by ExactSolution::_compute_error(), GMVIO::copy_nodal_solution(), ExactErrorEstimator::estimate_error(), ExactErrorEstimator::find_squared_element_error(), System::variable_type(), EnsightIO::write_scalar_ascii(), and EnsightIO::write_vector_ascii().

{
  // Make sure the variable exists
  std::map<std::string, unsigned short int>::const_iterator
    pos = _variable_numbers.find(var);
  
  if (pos == _variable_numbers.end())
    {
      std::cerr << 'ERROR: variable '
                << var
                << ' does not exist in this system!'
                << std::endl;      
      libmesh_error();
    }
  libmesh_assert (_variables[pos->second].name() == var);

  return pos->second;
}
 

const FEType & System::variable_type (const unsigned inti) const [inline, inherited]Returns:

the finite element type variable number i.

Definition at line 1147 of file system.h.

References System::_variables.

Referenced by System::add_variable(), EquationSystems::build_discontinuous_solution_vector(), EquationSystems::build_solution_vector(), GMVIO::copy_nodal_solution(), FEMContext::FEMContext(), System::write_header(), EnsightIO::write_scalar_ascii(), and EnsightIO::write_vector_ascii().

{
  libmesh_assert (i < _variables.size());
  
  return _variables[i].type();
}
 

const FEType & System::variable_type (const std::string &var) const [inline, inherited]Returns:

the finite element type for variable var.

Definition at line 1157 of file system.h.

References System::_variables, and System::variable_number().

{
  return _variables[this->variable_number(var)].type();
}
 

const std::string & System::vector_name (const unsigned intvec_num) [inherited]Returns:

the name of this system's additional vector number vec_num (where the vectors are counted starting with 0).

Definition at line 621 of file system.C.

References System::vectors_begin(), and System::vectors_end().

{
  const_vectors_iterator v = vectors_begin();
  const_vectors_iterator v_end = vectors_end();
  unsigned int num = 0;
  while((num<vec_num) && (v!=v_end))
    {
      num++;
      ++v;
    }
  libmesh_assert(v!=v_end);
  return v->first;
}
 

System::const_vectors_iterator System::vectors_begin () const [inline, inherited]Beginning of vectors container

Definition at line 1193 of file system.h.

References System::_vectors.

{
  return _vectors.begin();
}
 

System::vectors_iterator System::vectors_begin () [inline, inherited]Beginning of vectors container

Definition at line 1187 of file system.h.

References System::_vectors.

Referenced by UniformRefinementEstimator::_estimate_error(), System::get_vector(), VTKIO::system_vectors_to_vtk(), and System::vector_name().

{
  return _vectors.begin();
}
 

System::const_vectors_iterator System::vectors_end () const [inline, inherited]End of vectors container

Definition at line 1205 of file system.h.

References System::_vectors.

{
  return _vectors.end();
}
 

System::vectors_iterator System::vectors_end () [inline, inherited]End of vectors container

Definition at line 1199 of file system.h.

References System::_vectors.

Referenced by UniformRefinementEstimator::_estimate_error(), System::get_vector(), VTKIO::system_vectors_to_vtk(), and System::vector_name().

{
  return _vectors.end();
}
 

void System::write_header (Xdr &io, const std::string &version, const boolwrite_additional_data) const [inherited]Writes the basic data header for this System.

This method implements the output of a System object, embedded in the output of an EquationSystems<T_sys>. This warrants some documentation. The output of this part consists of 5 sections:

for this system

5.) The number of variables in the system (unsigned int)

for each variable in the system

6.) The name of the variable (string)

7.) Combined in an FEType:

The approximation order(s) of the variable (Order Enum, cast to int/s)
The finite element family/ies of the variable (FEFamily Enum, cast to int/s)

end variable loop

8.) The number of additional vectors (unsigned int),

for each additional vector in the system object

9.) the name of the additional vector (string)

end system

Definition at line 815 of file system_io.C.

References System::_vectors, Xdr::data(), FEType::family, System::get_mesh(), FEType::inf_map, System::n_vars(), System::n_vectors(), System::name(), FEType::order, MeshBase::processor_id(), FEType::radial_family, FEType::radial_order, System::variable_name(), System::variable_type(), and Xdr::writing().

{
  libmesh_assert (io.writing());


  // Only write the header information
  // if we are processor 0.
  if (this->get_mesh().processor_id() != 0)
    return;
  
  std::string comment;
  char buf[80];

  // 5.) 
  // Write the number of variables in the system

  {
    // set up the comment
    comment = '# No. of Variables in System '';
    comment += this->name();
    comment += ''';    
          
    unsigned int n_vars = this->n_vars();   
    io.data (n_vars, comment.c_str());
  }


  for (unsigned int var=0; var<this->n_vars(); var++)
    {
      // 6.)
      // Write the name of the var-th variable
      {
        // set up the comment
        comment  = '#   Name, Variable No. ';
        std::sprintf(buf, '%d', var);
        comment += buf;
        comment += ', System '';
        comment += this->name();
        comment += ''';
              
        std::string var_name = this->variable_name(var);             
        io.data (var_name, comment.c_str());
      }
      
      // 7.)
      // Write the approximation order of the var-th variable 
      // in this system
      {
        // set up the comment
        comment = '#     Approximation Order, Variable '';
        std::sprintf(buf, '%s', this->variable_name(var).c_str());
        comment += buf;
        comment += '', System '';
        comment += this->name();
        comment += ''';
              
        int order = static_cast<int>(this->variable_type(var).order);         
        io.data (order, comment.c_str());
      }


#ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
      
      // do the same for radial_order
      {
        comment = '#     Radial Approximation Order, Variable '';
        std::sprintf(buf, '%s', this->variable_name(var).c_str());
        comment += buf;
        comment += '', System '';
        comment += this->name();
        comment += ''';
      
        int rad_order = static_cast<int>(this->variable_type(var).radial_order);              
        io.data (rad_order, comment.c_str());
      }

#endif
      
      // Write the Finite Element type of the var-th variable 
      // in this System
      {
        // set up the comment
        comment = '#     FE Family, Variable '';
        std::sprintf(buf, '%s', this->variable_name(var).c_str());
        comment += buf;
        comment += '', System '';
        comment += this->name();
        comment += ''';
      
        const FEType& type = this->variable_type(var);        
        int fam = static_cast<int>(type.family);              
        io.data (fam, comment.c_str());

#ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
        
        comment = '#     Radial FE Family, Variable '';
        std::sprintf(buf, '%s', this->variable_name(var).c_str());
        comment += buf;
        comment += '', System '';
        comment += this->name();
        comment += ''';
      
        int radial_fam = static_cast<int>(type.radial_family);
        io.data (radial_fam, comment.c_str());  
        
        comment = '#     Infinite Mapping Type, Variable '';
        std::sprintf(buf, '%s', this->variable_name(var).c_str());
        comment += buf;
        comment += '', System '';
        comment += this->name();
        comment += ''';

        int i_map = static_cast<int>(type.inf_map);           
        io.data (i_map, comment.c_str());
#endif
      }
    } // end of the variable loop
 
  // 8.) 
  // Write the number of additional vectors in the System.
  // If write_additional_data==false, then write zero for
  // the number of additional vectors.
  {       
    {
      // set up the comment
      comment = '# No. of Additional Vectors, System '';
      comment += this->name();
      comment += ''';
    
      unsigned int n_vectors = write_additional_data ? this->n_vectors () : 0;
      io.data (n_vectors, comment.c_str());
    }  
      
    if (write_additional_data)
      {
        std::map<std::string, NumericVector<Number>* >::const_iterator
          vec_pos = this->_vectors.begin();
        unsigned int cnt=0;
        
        for (; vec_pos != this->_vectors.end(); ++vec_pos)
          {
            // 9.)
            // write the name of the cnt-th additional vector
            comment =  '# Name of ';
            std::sprintf(buf, '%d', cnt++);
            comment += buf;
            comment += 'th vector';
            std::string vec_name = vec_pos->first;
            
            io.data (vec_name, comment.c_str());
          }
      }
  }
}
 

void System::write_parallel_data (Xdr &io, const boolwrite_additional_data) const [inherited]Writes additional data, namely vectors, for this System. This method may safely be called on a distributed-memory mesh. This method will create an individual file for each processor in the simulation where the local solution components for that processor will be stored.

This method implements the output of the vectors contained in this System object, embedded in the output of an EquationSystems<T_sys>.

9.) The global solution vector, re-ordered to be node-major (More on this later.)

for each additional vector in the object

10.) The global additional vector, re-ordered to be node-major (More on this later.)

Note that the actual IO is handled through the Xdr class (to be renamed later?) which provides a uniform interface to both the XDR (eXternal Data Representation) interface and standard ASCII output. Thus this one section of code will read XDR or ASCII files with no changes.

This method implements the output of the vectors contained in this System object, embedded in the output of an EquationSystems<T_sys>.

9.) The global solution vector, re-ordered to be node-major (More on this later.)

for each additional vector in the object

10.) The global additional vector, re-ordered to be node-major (More on this later.)

Note that the actual IO is handled through the Xdr class (to be renamed later?) which provides a uniform interface to both the XDR (eXternal Data Representation) interface and standard ASCII output. Thus this one section of code will read XDR or ASCII files with no changes.

Definition at line 1213 of file system_io.C.

References System::_vectors, Xdr::data(), System::get_mesh(), DofObject::invalid_id, System::n_vars(), System::name(), System::number(), System::solution, and Xdr::writing().

{
  std::string comment;
  
  libmesh_assert (io.writing());

  std::vector<Number> io_buffer; io_buffer.reserve(this->solution->local_size());

  // build the ordered nodes and element maps.
  // when writing/reading parallel files we need to iterate
  // over our nodes/elements in order of increasing global id().
  // however, this is not guaranteed to be ordering we obtain
  // by using the node_iterators/element_iterators directly.
  // so build a set, sorted by id(), that provides the ordering.
  // further, for memory economy build the set but then transfer
  // its contents to vectors, which will be sorted.
  std::vector<const DofObject*> ordered_nodes, ordered_elements;
  {    
    std::set<const DofObject*, CompareDofObjectsByID>
      ordered_nodes_set (this->get_mesh().local_nodes_begin(),
                         this->get_mesh().local_nodes_end());
      
      ordered_nodes.insert(ordered_nodes.end(),
                           ordered_nodes_set.begin(),
                           ordered_nodes_set.end());
  }
  {
    std::set<const DofObject*, CompareDofObjectsByID>
      ordered_elements_set (this->get_mesh().local_elements_begin(),
                            this->get_mesh().local_elements_end());
      
      ordered_elements.insert(ordered_elements.end(),
                              ordered_elements_set.begin(),
                              ordered_elements_set.end());
  }
  
  const unsigned int sys_num = this->number();
  const unsigned int n_vars  = this->n_vars();
  
  // Loop over each variable and each node, and write out the value.
  for (unsigned int var=0; var<n_vars; var++)
    {
      // First write the node DOF values
      for (std::vector<const DofObject*>::const_iterator
             it = ordered_nodes.begin(); it != ordered_nodes.end(); ++it)      
        for (unsigned int comp=0; comp<(*it)->n_comp(sys_num, var); comp++)
          {
            //std::cout << '(*it)->id()=' << (*it)->id() << std::endl;
            libmesh_assert ((*it)->dof_number(sys_num, var, comp) !=
                    DofObject::invalid_id);
            
            io_buffer.push_back((*this->solution)((*it)->dof_number(sys_num, var, comp)));
          }

      // Then write the element DOF values
      for (std::vector<const DofObject*>::const_iterator
             it = ordered_elements.begin(); it != ordered_elements.end(); ++it)
        for (unsigned int comp=0; comp<(*it)->n_comp(sys_num, var); comp++)
          {
            libmesh_assert ((*it)->dof_number(sys_num, var, comp) !=
                    DofObject::invalid_id);
              
            io_buffer.push_back((*this->solution)((*it)->dof_number(sys_num, var, comp)));    
          }
    }
  
  // 9.)
  //
  // Actually write the reordered solution vector 
  // for the ith system to disk
  
  // set up the comment
  {
    comment = '# System '';
    comment += this->name();
    comment += '' Solution Vector';
  }
  
  io.data (io_buffer, comment.c_str());   
  
  // Only write additional vectors if wanted  
  if (write_additional_data)
    {     
      std::map<std::string, NumericVector<Number>* >::const_iterator
        pos = _vectors.begin();
  
      for(; pos != this->_vectors.end(); ++pos)
        {
          io_buffer.clear(); io_buffer.reserve( pos->second->local_size());
          
          // Loop over each variable and each node, and write out the value.
          for (unsigned int var=0; var<n_vars; var++)
            {
              // First write the node DOF values
              for (std::vector<const DofObject*>::const_iterator
                     it = ordered_nodes.begin(); it != ordered_nodes.end(); ++it)      
                for (unsigned int comp=0; comp<(*it)->n_comp(sys_num, var); comp++)
                  {
                    libmesh_assert ((*it)->dof_number(sys_num, var, comp) !=
                            DofObject::invalid_id);
                      
                    io_buffer.push_back((*pos->second)((*it)->dof_number(sys_num, var, comp)));   
                  }

              // Then write the element DOF values
              for (std::vector<const DofObject*>::const_iterator
                     it = ordered_elements.begin(); it != ordered_elements.end(); ++it)
                for (unsigned int comp=0; comp<(*it)->n_comp(sys_num, var); comp++)
                  {
                    libmesh_assert ((*it)->dof_number(sys_num, var, comp) !=
                            DofObject::invalid_id);
              
                    io_buffer.push_back((*pos->second)((*it)->dof_number(sys_num, var, comp)));
                  }
            }
          
          // 10.)
          //
          // Actually write the reordered additional vector 
          // for this system to disk
          
          // set up the comment
          {
            comment = '# System '';
            comment += this->name(); 
            comment += '' Additional Vector '';
            comment += pos->first;
            comment += ''';
          }
              
          io.data (io_buffer, comment.c_str());         
        }
    }
}
 

void System::write_serialized_data (Xdr &io, const boolwrite_additional_data = true) const [inherited]Writes additional data, namely vectors, for this System. This method may safely be called on a distributed-memory mesh.

This method implements the output of the vectors contained in this System object, embedded in the output of an EquationSystems<T_sys>.

9.) The global solution vector, re-ordered to be node-major (More on this later.)

for each additional vector in the object

10.) The global additional vector, re-ordered to be node-major (More on this later.)

Definition at line 1370 of file system_io.C.

References System::_vectors, Xdr::comment(), System::name(), libMesh::processor_id(), System::solution, and System::write_serialized_vector().

{
  parallel_only();
  std::string comment;

  this->write_serialized_vector(io, *this->solution); 

  // set up the comment
  if (libMesh::processor_id() == 0)
    {
      comment = '# System '';
      comment += this->name();
      comment += '' Solution Vector';

      io.comment (comment);
    }

  // Only write additional vectors if wanted  
  if (write_additional_data)
    {     
      std::map<std::string, NumericVector<Number>* >::const_iterator
        pos = _vectors.begin();
  
      for(; pos != this->_vectors.end(); ++pos)
        {
          this->write_serialized_vector(io, *pos->second);

          // set up the comment
          if (libMesh::processor_id() == 0)
            {
              comment = '# System '';
              comment += this->name(); 
              comment += '' Additional Vector '';
              comment += pos->first;
              comment += ''';
              io.comment (comment);       
            }
        }
    }
}
 

void System::zero_variable (NumericVector< Number > &v, unsigned intvar_num) const [inherited]Zeroes all dofs in v that correspond to variable number var_num.

Definition at line 763 of file system.C.

References MeshBase::active_local_elements_begin(), MeshBase::active_local_elements_end(), DofObject::dof_number(), System::get_mesh(), MeshBase::local_nodes_begin(), MeshBase::local_nodes_end(), DofObject::n_comp(), System::n_vars(), System::number(), and NumericVector< T >::set().

{
  /* Make sure the call makes sense.  */
  libmesh_assert(var_num<this->n_vars());

  /* Get a reference to the mesh.  */
  const MeshBase& mesh = this->get_mesh();

  /* Check which system we are.  */
  const unsigned int sys_num = this->number();

  /* Loop over nodes.  */
  {
    MeshBase::const_node_iterator it = mesh.local_nodes_begin();
    const MeshBase::const_node_iterator end_it = mesh.local_nodes_end(); 
    for ( ; it != end_it; ++it)
      {    
        const Node* node = *it;
        unsigned int n_comp = node->n_comp(sys_num,var_num);
        for(unsigned int i=0; i<n_comp; i++)
          {
            const unsigned int index = node->dof_number(sys_num,var_num,i);
            v.set(index,0.0);
          }
      }
  }

  /* Loop over elements.  */
  {
    MeshBase::const_element_iterator it = mesh.active_local_elements_begin();
    const MeshBase::const_element_iterator end_it = mesh.active_local_elements_end(); 
    for ( ; it != end_it; ++it)
      {    
        const Elem* elem = *it;
        unsigned int n_comp = elem->n_comp(sys_num,var_num);
        for(unsigned int i=0; i<n_comp; i++)
          {
            const unsigned int index = elem->dof_number(sys_num,var_num,i);
            v.set(index,0.0);
          }
      }
  }
}
 

Member Data Documentation

 

ReferenceCounter::Counts ReferenceCounter::_counts [static, protected, inherited]Actually holds the data.

Definition at line 110 of file reference_counter.h.

Referenced by ReferenceCounter::get_info(), ReferenceCounter::increment_constructor_count(), and ReferenceCounter::increment_destructor_count().  

unsigned int FEMSystem::_mesh_sys [protected, inherited]System and variables from which to acquire moving mesh information

Definition at line 307 of file fem_system.h.

Referenced by FEMSystem::eulerian_residual(), FEMSystem::mesh_position_get(), FEMSystem::mesh_position_set(), FEMSystem::mesh_x_position(), FEMSystem::mesh_y_position(), FEMSystem::mesh_z_position(), and FEMSystem::numerical_jacobian().  

unsigned int FEMSystem::_mesh_x_var [protected, inherited]

Definition at line 307 of file fem_system.h.

Referenced by FEMSystem::eulerian_residual(), FEMSystem::mesh_position_get(), FEMSystem::mesh_x_position(), and FEMSystem::numerical_jacobian().  

unsigned int FEMSystem::_mesh_y_var [protected, inherited]

Definition at line 307 of file fem_system.h.

Referenced by FEMSystem::eulerian_residual(), FEMSystem::mesh_position_get(), FEMSystem::mesh_y_position(), and FEMSystem::numerical_jacobian().  

unsigned int FEMSystem::_mesh_z_var [protected, inherited]

Definition at line 307 of file fem_system.h.

Referenced by FEMSystem::eulerian_residual(), FEMSystem::mesh_position_get(), FEMSystem::mesh_z_position(), and FEMSystem::numerical_jacobian().  

Threads::spin_mutex ReferenceCounter::_mutex [static, protected, inherited]Mutual exclusion object to enable thread-safe reference counting.

Definition at line 123 of file reference_counter.h.  

Threads::atomic< unsigned int > ReferenceCounter::_n_objects [static, protected, inherited]The number of objects. Print the reference count information when the number returns to 0.

Definition at line 118 of file reference_counter.h.

Referenced by ReferenceCounter::n_objects(), ReferenceCounter::ReferenceCounter(), and ReferenceCounter::~ReferenceCounter().  

std::vector<bool> DifferentiableSystem::_time_evolving [protected, inherited]Stores bools to tell us which variables are evolving in time and which are just constraints

Definition at line 424 of file diff_system.h.

Referenced by DifferentiableSystem::clear(), FEMSystem::eulerian_residual(), FEMSystem::init_context(), DifferentiableSystem::init_data(), FEMSystem::mass_residual(), and DifferentiableSystem::time_evolving().  

bool System::assemble_before_solve [inherited]Flag which tells the system to whether or not to call the user assembly function during each call to solve(). By default, every call to solve() begins with a call to the user assemble, so this flag is true. (For explicit systems, 'solving' the system occurs during the assembly step, so this flag is always true for explicit systems.)

You will only want to set this to false if you need direct control over when the system is assembled, and are willing to track the state of its assembly yourself. An example of such a case is an implicit system with multiple right hand sides. In this instance, a single assembly would likely be followed with multiple calls to solve.

The frequency system and Newmark system have their own versions of this flag, called _finished_assemble, which might be able to be replaced with this more general concept.

Definition at line 713 of file system.h.

Referenced by NonlinearImplicitSystem::adjoint_solve(), LinearImplicitSystem::adjoint_solve(), LinearImplicitSystem::solve(), and EigenSystem::solve().  

bool DifferentiableSystem::assemble_qoi_sides [inherited]If assemble_qoi_sides is true (it is false by default), the assembly loop for a quantity of interest or its derivatives will loop over all sides as well as all elements.

Definition at line 199 of file diff_system.h.

Referenced by FEMSystem::assemble_qoi(), and FEMSystem::assemble_qoi_derivative().  

bool DifferentiableSystem::compute_internal_sides [inherited]compute_internal_sides is false by default, indicating that side_* computations will only be done on boundary sides. If compute_internal_sides is true, computations will be done on sides between elements as well.

Definition at line 147 of file diff_system.h.

Referenced by FEMSystem::assemble_qoi(), FEMSystem::assemble_qoi_derivative(), FEMSystem::assembly(), and FEMSystem::postprocess().  

Real* ContinuationSystem::continuation_parameterThe continuation parameter must be a member variable of the derived class, and the 'continuation_parameter' pointer defined here must be a pointer to that variable. This is how the continuation system updates the derived class's continuation parameter.

Also sometimes referred to as 'lambda' in the code comments.

Definition at line 115 of file continuation_system.h.

Referenced by apply_predictor(), continuation_solve(), initialize_tangent(), save_current_solution(), solve_tangent(), and update_solution().  

Real ContinuationSystem::continuation_parameter_toleranceHow tightly should the Newton iterations attempt to converge delta_lambda. Defaults to 1.e-6.

Definition at line 134 of file continuation_system.h.

Referenced by continuation_solve().  

AutoPtr<NumericVector<Number> > System::current_local_solution [inherited]All the values I need to compute my contribution to the simulation at hand. Think of this as the current solution with any ghost values needed from other processors. This vector is necessarily larger than the solution vector in the case of a parallel simulation. The update() member is used to synchronize the contents of the solution and current_local_solution vectors.

Definition at line 740 of file system.h.

Referenced by UniformRefinementEstimator::_estimate_error(), NonlinearImplicitSystem::adjoint_solve(), System::clear(), Problem_Interface::computeF(), System::current_solution(), ExactErrorEstimator::estimate_error(), System::init_data(), System::project_solution(), System::re_update(), System::reinit(), System::restrict_vectors(), and System::update().  

NumericVector<Number>* ContinuationSystem::delta_u [private]Temporary vector 'delta u' ... the Newton step update in our custom augmented PDE solve.

Definition at line 370 of file continuation_system.h.

Referenced by continuation_solve(), init_data(), solve_tangent(), and update_solution().  

Real DifferentiableSystem::deltat [inherited]For time-dependent problems, this is the amount delta t to advance the solution in time.

Definition at line 365 of file diff_system.h.

Referenced by UnsteadySolver::advance_timestep(), EulerSolver::element_residual(), Euler2Solver::element_residual(), FEMSystem::eulerian_residual(), EulerSolver::side_residual(), Euler2Solver::side_residual(), UnsteadySolver::solve(), and TwostepTimeSolver::solve().  

Real ContinuationSystem::dlambda_ds [private]The most recent value of the derivative of the continuation parameter with respect to s. We use 'lambda' here for shortness of notation, lambda always refers to the continuation parameter.

Definition at line 396 of file continuation_system.h.

Referenced by apply_predictor(), continuation_solve(), initialize_tangent(), set_Theta_LOCA(), and solve_tangent().  

Real ContinuationSystem::ds [private]The initial arclength stepsize, selected by the user. This is the max-allowable arclength stepsize, but the algorithm may frequently reduce ds near turning points.

Definition at line 403 of file continuation_system.h.

Referenced by set_max_arclength_stepsize(), and update_solution().  

Real ContinuationSystem::ds_current [private]Value of stepsize currently in use. Will not exceed user-provided maximum arclength stepize ds.

Definition at line 409 of file continuation_system.h.

Referenced by apply_predictor(), continuation_solve(), initialize_tangent(), set_max_arclength_stepsize(), and update_solution().  

Real ContinuationSystem::ds_minThe minimum-allowed steplength, defaults to 1.e-8.

Definition at line 214 of file continuation_system.h.

Referenced by update_solution().  

NumericVector<Number>* ContinuationSystem::du_ds [private]Extra work vectors used by the continuation algorithm. These are added to the system by the init_data() routine.

The 'solution' tangent vector du/ds.

Definition at line 338 of file continuation_system.h.

Referenced by apply_predictor(), continuation_solve(), init_data(), initialize_tangent(), and solve_tangent().  

int FEMSystem::extra_quadrature_order [inherited]By default, when calling the user-defined residual functions, the FEMSystem will first set up an appropriate FEType::default_quadrature_rule() object for performing the integration. This rule will integrate elements of order up to 2*p+1 exactly (where p is the sum of the base FEType and local p refinement levels), but if additional (or reduced) quadrature accuracy is desired then this extra_quadrature_order (default 0) will be added.

Definition at line 247 of file fem_system.h.

Referenced by FEMContext::FEMContext().  

bool FEMSystem::fe_reinit_during_postprocess [inherited]If fe_reinit_during_postprocess is true (it is true by default), FE objects will be reinit()ed with their default quadrature rules. If false, FE objects will need to be reinit()ed by the user or will be in an undefined state.

Definition at line 236 of file fem_system.h.

Referenced by FEMSystem::postprocess().  

Real ContinuationSystem::initial_newton_toleranceHow much to try to reduce the residual by at the first (inexact) Newton step. This is frequently something large like 1/2 in an inexact Newton method, to prevent oversolving.

Definition at line 147 of file continuation_system.h.

Referenced by continuation_solve().  

AutoPtr<LinearSolver<Number> > ContinuationSystem::linear_solver [private]We maintain our own linear solver interface, for solving custom systems of equations and/or things which do not require a full-blown NewtonSolver.

Definition at line 377 of file continuation_system.h.

Referenced by continuation_solve(), and solve_tangent().  

SparseMatrix<Number>* ImplicitSystem::matrix [inherited]The system matrix. Implicit systems are characterized by the need to solve the linear system Ax=b. This is the system matrix A.

Definition at line 156 of file implicit_system.h.

Referenced by __libmesh_petsc_diff_solver_jacobian(), ImplicitSystem::add_system_matrix(), PetscDiffSolver::adjoint_solve(), NonlinearImplicitSystem::adjoint_solve(), NewtonSolver::adjoint_solve(), LinearImplicitSystem::adjoint_solve(), ImplicitSystem::assemble(), FEMSystem::assembly(), ImplicitSystem::clear(), NewmarkSystem::compute_matrix(), continuation_solve(), ImplicitSystem::init_matrices(), PetscDiffSolver::solve(), NonlinearImplicitSystem::solve(), NewtonSolver::solve(), LinearImplicitSystem::solve(), FrequencySystem::solve(), EigenTimeSolver::solve(), and solve_tangent().  

Real ContinuationSystem::max_continuation_parameterThe maximum-allowable value of the continuation parameter. The Newton iterations will quit if the continuation parameter goes above this value. If this value is zero, there is no maximum value for the continuation parameter.

Definition at line 174 of file continuation_system.h.

Referenced by continuation_solve().  

Real ContinuationSystem::min_continuation_parameterThe minimum-allowable value of the continuation parameter. The Newton iterations will quit if the continuation parameter falls below this value.

Definition at line 167 of file continuation_system.h.

Referenced by continuation_solve().  

unsigned int ContinuationSystem::n_arclength_reductionsNumber of arclength reductions to try when Newton fails to reduce the residual. For each arclength reduction, the arcstep size is cut in half.

Definition at line 209 of file continuation_system.h.

Referenced by continuation_solve().  

unsigned int ContinuationSystem::n_backtrack_stepsAnother scaling parameter suggested by the LOCA people. This one attempts to shrink the stepsize ds whenever the angle between the previous two tangent vectors gets large. Number of (Newton) backtracking steps to attempt if a Newton step does not reduce the residual. This is backtracking within a *single* Newton step, if you want to try a smaller arcstep, set n_arclength_reductions > 0.

Definition at line 202 of file continuation_system.h.

Referenced by continuation_solve().  

bool ContinuationSystem::newton_progress_checkTrue by default, the Newton progress check allows the Newton loop to exit if half the allowed iterations have elapsed without a reduction in the *initial* residual. In our experience this usually means the Newton steps are going to fail eventually and we could save some time by quitting early.

Definition at line 256 of file continuation_system.h.

Referenced by continuation_solve().  

NewtonSolver* ContinuationSystem::newton_solver [private]A pointer to the underlying Newton solver used by the DiffSystem. From this pointer, we can get access to all the parameters and options which are available to the 'normal' Newton solver.

Definition at line 389 of file continuation_system.h.

Referenced by continuation_solve(), solve_tangent(), and update_solution().  

unsigned int ContinuationSystem::newton_step [private]Loop counter for nonlinear (Newton) iteration loop.

Definition at line 424 of file continuation_system.h.

Referenced by continuation_solve(), and update_solution().  

Real ContinuationSystem::newton_stepgrowth_aggressivenessA measure of how rapidly one should attempt to grow the arclength stepsize based on the number of Newton iterations required to solve the problem. Default value is 1.0, if set to zero, will not try to grow or shrink the arclength stepsize based on the number of Newton iterations required.

Definition at line 247 of file continuation_system.h.

Referenced by update_solution().  

Real FEMSystem::numerical_jacobian_h [inherited]If calculating numeric jacobians is required, the FEMSystem will perturb each solution vector entry by numerical_jacobian_h when calculating finite differences.

Definition at line 254 of file fem_system.h.

Referenced by FEMSystem::numerical_jacobian().  

Real ContinuationSystem::old_continuation_parameterThe system also keeps track of the old value of the continuation parameter.

Definition at line 161 of file continuation_system.h.

Referenced by continuation_solve(), initialize_tangent(), save_current_solution(), solve_tangent(), and update_solution().  

bool DifferentiableSystem::postprocess_sides [inherited]If postprocess_sides is true (it is false by default), the postprocessing loop will loop over all sides as well as all elements.

Definition at line 192 of file diff_system.h.

Referenced by FEMSystem::postprocess().  

Predictor ContinuationSystem::predictor

Definition at line 238 of file continuation_system.h.

Referenced by apply_predictor().  

Real ContinuationSystem::previous_dlambda_ds [private]The old parameter tangent value.

Definition at line 414 of file continuation_system.h.

Referenced by apply_predictor(), and solve_tangent().  

Real ContinuationSystem::previous_ds [private]The previous arcstep length used.

Definition at line 419 of file continuation_system.h.

Referenced by apply_predictor(), and update_solution().  

NumericVector<Number>* ContinuationSystem::previous_du_ds [private]The value of du_ds from the previous solution

Definition at line 343 of file continuation_system.h.

Referenced by apply_predictor(), init_data(), and solve_tangent().  

NumericVector<Number>* ContinuationSystem::previous_u [private]The solution at the previous value of the continuation variable.

Definition at line 348 of file continuation_system.h.

Referenced by continuation_solve(), init_data(), initialize_tangent(), save_current_solution(), solve_tangent(), and update_solution().  

bool DifferentiableSystem::print_element_jacobians [inherited]Set print_element_jacobians to true to print each J_elem contribution.

Definition at line 402 of file diff_system.h.

Referenced by FEMSystem::assembly().  

bool DifferentiableSystem::print_jacobian_norms [inherited]Set print_jacobian_norms to true to print |J| whenever it is assembled.

Definition at line 392 of file diff_system.h.

Referenced by FEMSystem::assembly().  

bool DifferentiableSystem::print_jacobians [inherited]Set print_jacobians to true to print J whenever it is assembled.

Definition at line 397 of file diff_system.h.

Referenced by FEMSystem::assembly().  

bool DifferentiableSystem::print_residual_norms [inherited]Set print_residual_norms to true to print |F| whenever it is assembled.

Definition at line 382 of file diff_system.h.

Referenced by FEMSystem::assembly().  

bool DifferentiableSystem::print_residuals [inherited]Set print_residuals to true to print F whenever it is assembled.

Definition at line 387 of file diff_system.h.

Referenced by FEMSystem::assembly().  

bool DifferentiableSystem::print_solution_norms [inherited]Set print_residual_norms to true to print |U| whenever it is used in an assembly() call

Definition at line 371 of file diff_system.h.

Referenced by FEMSystem::assembly().  

bool DifferentiableSystem::print_solutions [inherited]Set print_solutions to true to print U whenever it is used in an assembly() call

Definition at line 377 of file diff_system.h.

Referenced by FEMSystem::assembly().  

Number System::qoi [inherited]Value of the quantity of interest

Definition at line 745 of file system.h.

Referenced by FEMSystem::assemble_qoi(), and ExplicitSystem::assemble_qoi().  

bool ContinuationSystem::quietIf quiet==false, the System prints extra information about what it is doing.

Definition at line 121 of file continuation_system.h.

Referenced by continuation_solve(), initialize_tangent(), set_Theta(), set_Theta_LOCA(), solve_tangent(), and update_solution().  

NumericVector<Number>* ExplicitSystem::rhs [inherited]The system matrix. Implicit systems are characterized by the need to solve the linear system Ax=b. This is the right-hand-side vector b.

Definition at line 132 of file explicit_system.h.

Referenced by __libmesh_petsc_diff_solver_residual(), ExplicitSystem::add_system_rhs(), PetscDiffSolver::adjoint_solve(), NonlinearImplicitSystem::adjoint_solve(), NewtonSolver::adjoint_solve(), LinearImplicitSystem::adjoint_solve(), ImplicitSystem::assemble(), ExplicitSystem::assemble_qoi(), FEMSystem::assemble_qoi_derivative(), ExplicitSystem::assemble_qoi_derivative(), FEMSystem::assembly(), ExplicitSystem::clear(), continuation_solve(), NewtonSolver::line_search(), PetscDiffSolver::solve(), NonlinearImplicitSystem::solve(), NewtonSolver::solve(), LinearImplicitSystem::solve(), FrequencySystem::solve(), solve_tangent(), and NewmarkSystem::update_rhs().  

RHS_Mode ContinuationSystem::rhs_mode [protected]

Definition at line 289 of file continuation_system.h.

Referenced by continuation_solve(), solve(), and solve_tangent().  

AutoPtr<NumericVector<Number> > System::solution [inherited]Data structure to hold solution values.

Definition at line 728 of file system.h.

Referenced by __libmesh_petsc_diff_solver_jacobian(), __libmesh_petsc_diff_solver_residual(), ExactSolution::_compute_error(), UniformRefinementEstimator::_estimate_error(), PetscDiffSolver::adjoint_solve(), UnsteadySolver::advance_timestep(), AdaptiveTimeSolver::advance_timestep(), apply_predictor(), FEMSystem::assembly(), System::clear(), System::compare(), Problem_Interface::computeF(), continuation_solve(), GMVIO::copy_nodal_solution(), UnsteadySolver::du(), DofMap::enforce_constraints_exactly(), PatchRecoveryErrorEstimator::estimate_error(), JumpErrorEstimator::estimate_error(), ExactErrorEstimator::estimate_error(), EigenSystem::get_eigenpair(), System::init_data(), initialize_tangent(), DofMap::max_constraint_error(), FEMSystem::mesh_position_get(), ErrorVector::plot_error(), System::project_solution(), System::re_update(), System::read_legacy_data(), System::read_parallel_data(), System::read_serialized_data(), System::reinit(), System::restrict_vectors(), save_current_solution(), VTKIO::solution_to_vtk(), TwostepTimeSolver::solve(), PetscDiffSolver::solve(), NonlinearImplicitSystem::solve(), NewtonSolver::solve(), LinearImplicitSystem::solve(), FrequencySystem::solve(), solve_tangent(), System::update(), System::update_global_solution(), update_solution(), NewmarkSystem::update_u_v_a(), System::write_parallel_data(), and System::write_serialized_data().  

Real ContinuationSystem::solution_toleranceHow tightly should the Newton iterations attempt to converge ||delta_u|| Defaults to 1.e-6.

Definition at line 140 of file continuation_system.h.

Referenced by continuation_solve().  

bool ContinuationSystem::tangent_initialized [private]False until initialize_tangent() is called

Definition at line 382 of file continuation_system.h.

Referenced by continuation_solve(), initialize_tangent(), solve_tangent(), and update_solution().  

Real ContinuationSystem::ThetaArclength normalization parameter. Defaults to 1.0 (no normalization). Used to ensure that one term in the arclength contstraint equation does not wash out all the others.

Definition at line 181 of file continuation_system.h.

Referenced by continuation_solve(), initialize_tangent(), set_Theta(), solve_tangent(), and update_solution().  

Real ContinuationSystem::Theta_LOCAAnother normalization parameter, which is described in the LOCA manual. This one is designed to maintain a relatively 'fixed' value of d(lambda)/ds. It is initially 1.0 and is updated after each solve.

Definition at line 188 of file continuation_system.h.

Referenced by continuation_solve(), initialize_tangent(), set_Theta_LOCA(), solve_tangent(), and update_solution().  

Real DifferentiableSystem::time [inherited]For time-dependent problems, this is the time t at the beginning of the current timestep. But do *not* access this time during an assembly! Use the DiffContext::time value instead to get correct results.

Definition at line 359 of file diff_system.h.

Referenced by UnsteadySolver::advance_timestep(), AdaptiveTimeSolver::advance_timestep(), and TwostepTimeSolver::solve().  

AutoPtr<TimeSolver> DifferentiableSystem::time_solver [inherited]A pointer to the solver object we're going to use. This must be instantiated by the user before solving!

Definition at line 351 of file diff_system.h.

Referenced by DifferentiableSystem::adjoint_solve(), FEMSystem::assembly(), continuation_solve(), FEMSystem::eulerian_residual(), DifferentiableSystem::init_data(), DifferentiableSystem::reinit(), DifferentiableSystem::solve(), and solve_tangent().  

bool DifferentiableSystem::use_fixed_solution [inherited]A boolean to be set to true by systems using elem_fixed_solution, so that the library code will create it. False by default. Warning: if this variable is set to true, it must be before init_data() is called.

Definition at line 411 of file diff_system.h.

Referenced by DifferentiableSystem::clear(), DiffContext::DiffContext(), SteadySolver::element_residual(), EulerSolver::element_residual(), Euler2Solver::element_residual(), FEMContext::reinit(), SteadySolver::side_residual(), EulerSolver::side_residual(), and Euler2Solver::side_residual().  

Real FEMSystem::verify_analytic_jacobians [inherited]If verify_analytic_jacobian is equal to zero (as it is by default), no numeric jacobians will be calculated unless an overloaded element_time_derivative(), element_constraint(), side_time_derivative(), or side_constraint() function cannot provide an analytic jacobian upon request.

If verify_analytic_jacobian is equal to the positive value tol, then any time a full analytic element jacobian can be calculated it will be tested against a numerical jacobian on the same element, and the program will abort if the relative error (in matrix l1 norms) exceeds tol.

Definition at line 269 of file fem_system.h.

Referenced by FEMSystem::assembly().  

NumericVector<Number>* ContinuationSystem::y [private]Temporary vector 'y' ... stores -du/dlambda, the solution of Ay=G_{}.

Definition at line 353 of file continuation_system.h.

Referenced by continuation_solve(), init_data(), initialize_tangent(), solve_tangent(), and update_solution().  

NumericVector<Number>* ContinuationSystem::y_old [private]Temporary vector 'y_old' ... stores the previous value of -du/dlambda, which is the solution of Ay=G_{}.

Definition at line 359 of file continuation_system.h.

Referenced by continuation_solve(), init_data(), and update_solution().  

NumericVector<Number>* ContinuationSystem::z [private]Temporary vector 'z' ... the solution of Az = -G

Definition at line 364 of file continuation_system.h.

Referenced by continuation_solve(), and init_data().

 

Author

Generated automatically by Doxygen for libMesh from the source code.


 

Index

NAME
SYNOPSIS
Public Types
Public Member Functions
Static Public Member Functions
Public Attributes
Protected Types
Protected Member Functions
Protected Attributes
Static Protected Attributes
Private Member Functions
Private Attributes
Detailed Description
Member Typedef Documentation
typedef std::map<std::string, SparseMatrix<Number>* >::const_iterator ImplicitSystem::const_matrices_iterator [inherited]
typedef std::map<std::string, NumericVector<Number>* >::const_iterator System::const_vectors_iterator [inherited]
typedef std::map<std::string, std::pair<unsigned int, unsigned int> > ReferenceCounter::Counts [protected, inherited]Data structure to log the information. The log is identified by the class name.
typedef std::map<std::string, SparseMatrix<Number>* >::iterator ImplicitSystem::matrices_iterator [inherited]Matrix iterator typedefs.
typedef FEMSystem ContinuationSystem::ParentThe type of the parent.
typedef ContinuationSystem ContinuationSystem::sys_typeThe type of system.
typedef bool(TimeSolver::* FEMSystem::TimeSolverResPtr)(bool, DiffContext &) [protected, inherited]Syntax sugar to make numerical_jacobian() declaration easier.
typedef std::map<std::string, NumericVector<Number>* >::iterator System::vectors_iterator [inherited]Vector iterator typedefs.
Member Enumeration Documentation
enum ContinuationSystem::PredictorThe code provides the ability to select from different predictor schemes for getting the initial guess for the solution at the next point on the solution arc.
enum ContinuationSystem::RHS_Mode [protected]There are (so far) two different vectors which may be assembled using the assembly routine: 1.) The Residual = the normal PDE weighted residual 2.) G_Lambda = the derivative wrt the parameter lambda of the PDE weighted residual
Constructor & Destructor Documentation
ContinuationSystem::ContinuationSystem (EquationSystems &es, const std::string &name, const unsigned intnumber)Constructor. Optionally initializes required data structures.
ContinuationSystem::~ContinuationSystem () [virtual]Destructor.
Member Function Documentation
void System::activate () [inline, inherited]Activates the system. Only active systems are solved.
bool System::active () const [inline, inherited]Returns:
SparseMatrix< Number > & ImplicitSystem::add_matrix (const std::string &mat_name) [inherited]Adds the additional matrix mat_name to this system. Only allowed prior to assemble(). All additional matrices have the same sparsity pattern as the matrix used during solution. When not System but the user wants to initialize the mayor matrix, then all the additional matrices, if existent, have to be initialized by the user, too.
unsigned int System::add_variable (const std::string &var, const FEType &type, const std::set< subdomain_id_type > *constactive_subdomains = NULL) [inherited]Adds the variable var to the list of variables for this system. Returns the index number for the new variable.
unsigned int System::add_variable (const std::string &var, const Orderorder = FIRST, const FEFamilyfamily = LAGRANGE, const std::set< subdomain_id_type > *constactive_subdomains = NULL) [inherited]Adds the variable var to the list of variables for this system. Same as before, but assumes LAGRANGE as default value for FEType.family.
NumericVector< Number > & System::add_vector (const std::string &vec_name, const boolprojections = true) [inherited]Adds the additional vector vec_name to this system. Only allowed prior to init(). All the additional vectors are similarly distributed, like the solution, and inititialized to zero.
void DifferentiableSystem::adjoint_solve () [virtual, inherited]Invokes the adjoint solver associated with the system. For steady state solvers, this will find x where dF/dx = q. For transient solvers, this should integrate -dx/dt = dF/dx - q.
void ContinuationSystem::advance_arcstep ()Call this function after a continuation solve to compute the tangent and get the next initial guess.
void ContinuationSystem::apply_predictor () [private]Applies the predictor to the current solution to get a guess for the next solution.
void DifferentiableSystem::assemble () [virtual, inherited]Prepares matrix and rhs for matrix assembly. Users should not reimplement this
void FEMSystem::assemble_qoi () [virtual, inherited]Runs a qoi assembly loop over all elements, and if assemble_qoi_sides is true over all sides.
void FEMSystem::assemble_qoi_derivative () [virtual, inherited]Runs a qoi derivative assembly loop over all elements, and if assemble_qoi_sides is true over all sides.
void FEMSystem::assembly (boolget_residual, boolget_jacobian) [virtual, inherited]Prepares matrix or rhs for matrix assembly. Users may reimplement this to add pre- or post-assembly code before or after calling FEMSystem::assembly()
void System::attach_assemble_function (void fptrEquationSystems &es,const std::string &name) [inherited]Register a user function to use in assembling the system matrix and RHS.
void System::attach_constraint_function (void fptrEquationSystems &es,const std::string &name) [inherited]Register a user function for imposing constraints.
void System::attach_init_function (void fptrEquationSystems &es,const std::string &name) [inherited]Register a user function to use in initializing the system.
void System::attach_QOI_derivative (void fptrEquationSystems &es,const std::string &name) [inherited]Register a user function for evaluating derivatives of a quantity of interest with respect to test functions, whose values should be placed in System::rhs
void System::attach_QOI_function (void fptrEquationSystems &es,const std::string &name) [inherited]Register a user function for evaluating a quantity of interest, whose value should be placed in System::qoi
AutoPtr< DiffContext > FEMSystem::build_context () [virtual, inherited]Builds a FEMContext object with enough information to do evaluations on each element.
Real System::calculate_norm (NumericVector< Number > &v, unsigned intvar = 0, FEMNormTypenorm_type = L2) const [inherited]Returns:
Real System::calculate_norm (NumericVector< Number > &v, const SystemNorm &norm) const [inherited]Returns:
void ContinuationSystem::clear () [virtual]Clear all the data structures associated with the system.
bool System::compare (const System &other_system, const Realthreshold, const boolverbose) const [virtual, inherited]Returns:
void ContinuationSystem::continuation_solve ()Perform a continuation solve of the system. In general, you can only begin the continuation solves after either reading in or solving for two previous values of the control parameter. The prior two solutions are required for starting up the continuation method.
Number System::current_solution (const unsigned intglobal_dof_number) const [inherited]Returns:
void System::deactivate () [inline, inherited]Deactivates the system. Only active systems are solved.
virtual bool DifferentiableSystem::element_constraint (boolrequest_jacobian, DiffContext &) [inline, virtual, inherited]Adds the constraint contribution on elem to elem_residual. If this method receives request_jacobian = true, then it should compute elem_jacobian and return true if possible. If elem_jacobian has not been computed then the method should return false.
virtual void DifferentiableSystem::element_postprocess (DiffContext &) [inline, virtual, inherited]Does any work that needs to be done on elem in a postprocessing loop.
virtual void DifferentiableSystem::element_qoi (DiffContext &) [inline, virtual, inherited]Does any work that needs to be done on elem in a quantity of interest assembly loop, outputting to element_qoi.
virtual void DifferentiableSystem::element_qoi_derivative (DiffContext &) [inline, virtual, inherited]Does any work that needs to be done on elem in a quantity of interest derivative assembly loop, outputting to element_residual.
virtual bool DifferentiableSystem::element_time_derivative (boolrequest_jacobian, DiffContext &) [inline, virtual, inherited]Adds the time derivative contribution on elem to elem_residual. If this method receives request_jacobian = true, then it should compute elem_jacobian and return true if possible. If elem_jacobian has not been computed then the method should return false.
bool FEMSystem::eulerian_residual (boolrequest_jacobian, DiffContext &context) [virtual, inherited]Adds a pseudo-convection contribution on elem to elem_residual, if the nodes of elem are being translated by a moving mesh.
NumericVector< Number > & System::get_adjoint_solution () [inherited]Returns:
const NumericVector< Number > & System::get_adjoint_solution () const [inherited]Returns:
const DofMap & System::get_dof_map () const [inline, inherited]Returns:
DofMap & System::get_dof_map () [inline, inherited]Returns:
const EquationSystems& System::get_equation_systems () const [inline, inherited]Returns:
EquationSystems& System::get_equation_systems () [inline, inherited]Returns:
std::string System::get_info () const [inherited]Returns:
std::string ReferenceCounter::get_info () [static, inherited]Gets a string containing the reference information.
SparseMatrix< Number > & ImplicitSystem::get_matrix (const std::string &mat_name) [inherited]Returns:
const SparseMatrix< Number > & ImplicitSystem::get_matrix (const std::string &mat_name) const [inherited]Returns:
const MeshBase & System::get_mesh () const [inline, inherited]Returns:
MeshBase & System::get_mesh () [inline, inherited]Returns:
const NumericVector< Number > & System::get_vector (const std::string &vec_name) const [inherited]Returns:
NumericVector< Number > & System::get_vector (const std::string &vec_name) [inherited]Returns:
const NumericVector< Number > & System::get_vector (const unsigned intvec_num) const [inherited]Returns:
NumericVector< Number > & System::get_vector (const unsigned intvec_num) [inherited]Returns:
bool System::has_variable (const std::string &var) const [inherited]Returns:
bool ImplicitSystem::have_matrix (const std::string &mat_name) const [inline, inherited]Returns:
bool System::have_vector (const std::string &vec_name) const [inline, inherited]Returns:
void ReferenceCounter::increment_constructor_count (const std::string &name) [inline, protected, inherited]Increments the construction counter. Should be called in the constructor of any derived class that will be reference counted.
void ReferenceCounter::increment_destructor_count (const std::string &name) [inline, protected, inherited]Increments the destruction counter. Should be called in the destructor of any derived class that will be reference counted.
void System::init () [inherited]Initializes degrees of freedom on the current mesh. Sets the
void FEMSystem::init_context (DiffContext &c) [virtual, inherited]
void ContinuationSystem::init_data () [protected, virtual]Initializes the member data fields associated with the system, so that, e.g., assemble() may be used.
void ImplicitSystem::init_matrices () [protected, virtual, inherited]Initializes the matrices associated with this system.
void ContinuationSystem::initialize_tangent () [private]Before starting arclength continuation, we need at least 2 prior solutions (both solution and u_previous should be filled with meaningful values) And we need to initialize the tangent vector. This only needs to be called once.
void System::local_dof_indices (const unsigned intvar, std::set< unsigned int > &var_indices) const [inherited]Fills the std::set with the degrees of freedom on the local processor corresponding the the variable number passed in.
bool FEMSystem::mass_residual (boolrequest_jacobian, DiffContext &context) [virtual, inherited]Adds a mass vector contribution on elem to elem_residual. If this method receives request_jacobian = true, then it should compute elem_jacobian and return true if possible. If elem_jacobian has not been computed then the method should return false.
void FEMSystem::mesh_position_get () [inherited]Tells the FEMSystem to set the degree of freedom coefficients which should correspond to mesh nodal coordinates.
void FEMSystem::mesh_position_set () [inherited]Tells the FEMSystem to set the mesh nodal coordinates which should correspond to degree of freedom coefficients.
void FEMSystem::mesh_x_position (unsigned intsysnum, unsigned intvar) [virtual, inherited]Tells the FEMSystem that variable var from system number sysnum should be used to update the x coordinate of mesh nodes, in problems where the mesh itself is expected to move in time.
void FEMSystem::mesh_y_position (unsigned intsysnum, unsigned intvar) [virtual, inherited]Tells the FEMSystem that variable var from system number sysnum should be used to update the y coordinate of mesh nodes.
void FEMSystem::mesh_z_position (unsigned intsysnum, unsigned intvar) [virtual, inherited]Tells the FEMSystem that variable var from system number sysnum should be used to update the z coordinate of mesh nodes.
unsigned int System::n_active_dofs () const [inline, inherited]Returns the number of active degrees of freedom for this System.
unsigned int System::n_constrained_dofs () const [inherited]Returns:
unsigned int System::n_dofs () const [inherited]Returns:
unsigned int System::n_local_dofs () const [inherited]Returns:
unsigned int ImplicitSystem::n_matrices () const [inline, inherited]Returns:
static unsigned int ReferenceCounter::n_objects () [inline, static, inherited]Prints the number of outstanding (created, but not yet destroyed) objects.
unsigned int System::n_vars () const [inline, inherited]Returns:
unsigned int System::n_vectors () const [inline, inherited]Returns:
const std::string & System::name () const [inline, inherited]Returns:
unsigned int System::number () const [inline, inherited]Returns:
void FEMSystem::numerical_elem_jacobian (FEMContext &context) [protected, inherited]Uses the results of multiple element_residual() calls to numerically differentiate the corresponding jacobian on an element.
void FEMSystem::numerical_jacobian (TimeSolverResPtrres, FEMContext &context) [protected, inherited]Uses the results of multiple res calls to numerically differentiate the corresponding jacobian.
void FEMSystem::numerical_side_jacobian (FEMContext &context) [protected, inherited]Uses the results of multiple side_residual() calls to numerically differentiate the corresponding jacobian on an element's side.
void FEMSystem::postprocess () [virtual, inherited]Runs a postprocessing loop over all elements, and if postprocess_sides is true over all sides.
void ReferenceCounter::print_info () [static, inherited]Prints the reference information to std::cout.
void System::project_solution (Number fptrconst Point &p,const Parameters &parameters,const std::string &sys_name,const std::string &unknown_name, Gradient gptrconst Point &p,const Parameters &parameters,const std::string &sys_name,const std::string &unknown_name, Parameters &parameters) const [inherited]Projects the continuous functions onto the current solution.
bool& System::project_solution_on_reinit (void) [inline, inherited]Tells the System whether or not to project the solution vector onto new grids when the system is reinitialized. The solution will be projected unless project_solution_on_reinit() = false is called.
void System::project_vector (Number fptrconst Point &p,const Parameters &parameters,const std::string &sys_name,const std::string &unknown_name, Gradient gptrconst Point &p,const Parameters &parameters,const std::string &sys_name,const std::string &unknown_name, Parameters &parameters, NumericVector< Number > &new_vector) const [inherited]Projects the continuous functions onto the current mesh.
void System::project_vector (NumericVector< Number > &vector) const [protected, inherited]Projects the vector defined on the old mesh onto the new mesh.
void System::project_vector (const NumericVector< Number > &old_v, NumericVector< Number > &new_v) const [protected, inherited]Projects the vector defined on the old mesh onto the new mesh. The original vector is unchanged and the new vector is passed through the second argument.
void System::prolong_vectors () [virtual, inherited]Prolong vectors after the mesh has refined
void DifferentiableSystem::qoi_parameter_sensitivity (std::vector< Number * > &parameters, std::vector< Number > &sensitivities) [virtual, inherited]Solves for the derivative of the system's quantity of interest q with respect to each parameter p in parameters. Currently uses adjoint_solve, along with finite differenced derivatives (partial q / partial p) and (partial R / partial p).
void System::re_update () [virtual, inherited]Re-update the local values when the mesh has changed. This method takes the data updated by update() and makes it up-to-date on the current mesh.
void System::read_header (Xdr &io, const std::string &version, const boolread_header = true, const boolread_additional_data = true, const boolread_legacy_format = false) [inherited]Reads the basic data header for this System.
void System::read_legacy_data (Xdr &io, const boolread_additional_data = true) [inherited]Reads additional data, namely vectors, for this System.
void System::read_parallel_data (Xdr &io, const boolread_additional_data) [inherited]Reads additional data, namely vectors, for this System. This method may safely be called on a distributed-memory mesh. This method will read an individual file for each processor in the simulation where the local solution components for that processor are stored.
void System::read_serialized_data (Xdr &io, const boolread_additional_data = true) [inherited]Reads additional data, namely vectors, for this System. This method may safely be called on a distributed-memory mesh.
void DifferentiableSystem::reinit () [virtual, inherited]Reinitializes the member data fields associated with the system, so that, e.g., assemble() may be used.
void System::restrict_vectors () [virtual, inherited]Restrict vectors after the mesh has coarsened
void ContinuationSystem::save_current_solution ()Stores the current solution and continuation parameter (as 'previous_u' and 'old_continuation_paramter') for later referral. You may need to call this e.g. after the first regular solve, in order to store the first solution, before computing a second solution and beginning arclength continuation.
void ContinuationSystem::set_max_arclength_stepsize (Realmaxds) [inline]Sets (initializes) the max-allowable ds value and the current ds value. Call this before beginning arclength continuation. The default max stepsize is 0.1
void ContinuationSystem::set_Theta () [private]A centralized function for setting the normalization parameter theta
void ContinuationSystem::set_Theta_LOCA () [private]A centralized function for setting the other normalization parameter, i.e. the one suggested by the LOCA developers.
virtual bool DifferentiableSystem::side_constraint (boolrequest_jacobian, DiffContext &) [inline, virtual, inherited]Adds the time derivative contribution on side of elem to elem_residual. If this method receives request_jacobian = true, then it should compute elem_jacobian and return true if possible. If elem_jacobian has not been computed then the method should return false.
virtual bool DifferentiableSystem::side_mass_residual (boolrequest_jacobian, DiffContext &) [inline, virtual, inherited]Adds a mass vector contribution on side of elem to elem_residual. If this method receives request_jacobian = true, then it should compute elem_jacobian and return true if possible. If elem_jacobian has not been computed then the method should return false.
virtual void DifferentiableSystem::side_postprocess (DiffContext &) [inline, virtual, inherited]Does any work that needs to be done on side of elem in a postprocessing loop.
virtual void DifferentiableSystem::side_qoi (DiffContext &) [inline, virtual, inherited]Does any work that needs to be done on side of elem in a quantity of interest assembly loop, outputting to element_qoi.
virtual void DifferentiableSystem::side_qoi_derivative (DiffContext &) [inline, virtual, inherited]Does any work that needs to be done on side of elem in a quantity of interest derivative assembly loop, outputting to element_residual.
virtual bool DifferentiableSystem::side_time_derivative (boolrequest_jacobian, DiffContext &) [inline, virtual, inherited]Adds the time derivative contribution on side of elem to elem_residual. If this method receives request_jacobian = true, then it should compute elem_jacobian and return true if possible. If elem_jacobian has not been computed then the method should return false.
void ContinuationSystem::solve () [virtual]Perform a standard 'solve' of the system, without doing continuation.
void ContinuationSystem::solve_tangent () [private]Special solve algorithm for solving the tangent system.
sys_type& ImplicitSystem::system () [inline, inherited]Returns:
virtual std::string ImplicitSystem::system_type () const [inline, virtual, inherited]Assembles & solves the linear system Ax=b.
void FEMSystem::time_evolving (unsigned intvar) [virtual, inherited]Tells the FEMSystem that variable var is evolving with respect to time. In general, the user's init() function should call time_evolving() for any variables which behave like du/dt = F(u), and should not call time_evolving() for any variables which behave like 0 = G(u).
void System::update () [virtual, inherited]Update the local values to reflect the solution on neighboring processors.
void System::update_global_solution (std::vector< Number > &global_soln) const [inherited]Fill the input vector global_soln so that it contains the global solution on all processors. Requires communication with all other processors.
void System::update_global_solution (std::vector< Number > &global_soln, const unsigned intdest_proc) const [inherited]Fill the input vector global_soln so that it contains the global solution on processor dest_proc. Requires communication with all other processors.
void ContinuationSystem::update_solution () [private]This function (which assumes the most recent tangent vector has been computed) updates the solution and the control parameter with the initial guess for the next point on the continuation path.
void System::user_assembly () [virtual, inherited]Calls user's attached assembly function, or is overloaded by the user in derived classes.
void System::user_constrain () [virtual, inherited]Calls user's attached constraint function, or is overloaded by the user in derived classes.
void System::user_initialization () [virtual, inherited]Calls user's attached initialization function, or is overloaded by the user in derived classes.
void System::user_QOI () [virtual, inherited]Calls user's attached quantity of interest function, or is overloaded by the user in derived classes.
void System::user_QOI_derivative () [virtual, inherited]Calls user's attached quantity of interest derivative function, or is overloaded by the user in derived classes.
const System::Variable & System::variable (unsigned intvar) const [inline, inherited]Return a constant reference to Variable var.
const std::string & System::variable_name (const unsigned inti) const [inline, inherited]Returns:
unsigned short int System::variable_number (const std::string &var) const [inherited]Returns:
const FEType & System::variable_type (const unsigned inti) const [inline, inherited]Returns:
const FEType & System::variable_type (const std::string &var) const [inline, inherited]Returns:
const std::string & System::vector_name (const unsigned intvec_num) [inherited]Returns:
System::const_vectors_iterator System::vectors_begin () const [inline, inherited]Beginning of vectors container
System::vectors_iterator System::vectors_begin () [inline, inherited]Beginning of vectors container
System::const_vectors_iterator System::vectors_end () const [inline, inherited]End of vectors container
System::vectors_iterator System::vectors_end () [inline, inherited]End of vectors container
void System::write_header (Xdr &io, const std::string &version, const boolwrite_additional_data) const [inherited]Writes the basic data header for this System.
void System::write_parallel_data (Xdr &io, const boolwrite_additional_data) const [inherited]Writes additional data, namely vectors, for this System. This method may safely be called on a distributed-memory mesh. This method will create an individual file for each processor in the simulation where the local solution components for that processor will be stored.
void System::write_serialized_data (Xdr &io, const boolwrite_additional_data = true) const [inherited]Writes additional data, namely vectors, for this System. This method may safely be called on a distributed-memory mesh.
void System::zero_variable (NumericVector< Number > &v, unsigned intvar_num) const [inherited]Zeroes all dofs in v that correspond to variable number var_num.
Member Data Documentation
ReferenceCounter::Counts ReferenceCounter::_counts [static, protected, inherited]Actually holds the data.
unsigned int FEMSystem::_mesh_sys [protected, inherited]System and variables from which to acquire moving mesh information
unsigned int FEMSystem::_mesh_x_var [protected, inherited]
unsigned int FEMSystem::_mesh_y_var [protected, inherited]
unsigned int FEMSystem::_mesh_z_var [protected, inherited]
Threads::spin_mutex ReferenceCounter::_mutex [static, protected, inherited]Mutual exclusion object to enable thread-safe reference counting.
Threads::atomic< unsigned int > ReferenceCounter::_n_objects [static, protected, inherited]The number of objects. Print the reference count information when the number returns to 0.
std::vector<bool> DifferentiableSystem::_time_evolving [protected, inherited]Stores bools to tell us which variables are evolving in time and which are just constraints
bool System::assemble_before_solve [inherited]Flag which tells the system to whether or not to call the user assembly function during each call to solve(). By default, every call to solve() begins with a call to the user assemble, so this flag is true. (For explicit systems, 'solving' the system occurs during the assembly step, so this flag is always true for explicit systems.)
bool DifferentiableSystem::assemble_qoi_sides [inherited]If assemble_qoi_sides is true (it is false by default), the assembly loop for a quantity of interest or its derivatives will loop over all sides as well as all elements.
bool DifferentiableSystem::compute_internal_sides [inherited]compute_internal_sides is false by default, indicating that side_* computations will only be done on boundary sides. If compute_internal_sides is true, computations will be done on sides between elements as well.
Real* ContinuationSystem::continuation_parameterThe continuation parameter must be a member variable of the derived class, and the 'continuation_parameter' pointer defined here must be a pointer to that variable. This is how the continuation system updates the derived class's continuation parameter.
Real ContinuationSystem::continuation_parameter_toleranceHow tightly should the Newton iterations attempt to converge delta_lambda. Defaults to 1.e-6.
AutoPtr<NumericVector<Number> > System::current_local_solution [inherited]All the values I need to compute my contribution to the simulation at hand. Think of this as the current solution with any ghost values needed from other processors. This vector is necessarily larger than the solution vector in the case of a parallel simulation. The update() member is used to synchronize the contents of the solution and current_local_solution vectors.
NumericVector<Number>* ContinuationSystem::delta_u [private]Temporary vector 'delta u' ... the Newton step update in our custom augmented PDE solve.
Real DifferentiableSystem::deltat [inherited]For time-dependent problems, this is the amount delta t to advance the solution in time.
Real ContinuationSystem::dlambda_ds [private]The most recent value of the derivative of the continuation parameter with respect to s. We use 'lambda' here for shortness of notation, lambda always refers to the continuation parameter.
Real ContinuationSystem::ds [private]The initial arclength stepsize, selected by the user. This is the max-allowable arclength stepsize, but the algorithm may frequently reduce ds near turning points.
Real ContinuationSystem::ds_current [private]Value of stepsize currently in use. Will not exceed user-provided maximum arclength stepize ds.
Real ContinuationSystem::ds_minThe minimum-allowed steplength, defaults to 1.e-8.
NumericVector<Number>* ContinuationSystem::du_ds [private]Extra work vectors used by the continuation algorithm. These are added to the system by the init_data() routine.
int FEMSystem::extra_quadrature_order [inherited]By default, when calling the user-defined residual functions, the FEMSystem will first set up an appropriate FEType::default_quadrature_rule() object for performing the integration. This rule will integrate elements of order up to 2*p+1 exactly (where p is the sum of the base FEType and local p refinement levels), but if additional (or reduced) quadrature accuracy is desired then this extra_quadrature_order (default 0) will be added.
bool FEMSystem::fe_reinit_during_postprocess [inherited]If fe_reinit_during_postprocess is true (it is true by default), FE objects will be reinit()ed with their default quadrature rules. If false, FE objects will need to be reinit()ed by the user or will be in an undefined state.
Real ContinuationSystem::initial_newton_toleranceHow much to try to reduce the residual by at the first (inexact) Newton step. This is frequently something large like 1/2 in an inexact Newton method, to prevent oversolving.
AutoPtr<LinearSolver<Number> > ContinuationSystem::linear_solver [private]We maintain our own linear solver interface, for solving custom systems of equations and/or things which do not require a full-blown NewtonSolver.
SparseMatrix<Number>* ImplicitSystem::matrix [inherited]The system matrix. Implicit systems are characterized by the need to solve the linear system Ax=b. This is the system matrix A.
Real ContinuationSystem::max_continuation_parameterThe maximum-allowable value of the continuation parameter. The Newton iterations will quit if the continuation parameter goes above this value. If this value is zero, there is no maximum value for the continuation parameter.
Real ContinuationSystem::min_continuation_parameterThe minimum-allowable value of the continuation parameter. The Newton iterations will quit if the continuation parameter falls below this value.
unsigned int ContinuationSystem::n_arclength_reductionsNumber of arclength reductions to try when Newton fails to reduce the residual. For each arclength reduction, the arcstep size is cut in half.
unsigned int ContinuationSystem::n_backtrack_stepsAnother scaling parameter suggested by the LOCA people. This one attempts to shrink the stepsize ds whenever the angle between the previous two tangent vectors gets large. Number of (Newton) backtracking steps to attempt if a Newton step does not reduce the residual. This is backtracking within a *single* Newton step, if you want to try a smaller arcstep, set n_arclength_reductions > 0.
bool ContinuationSystem::newton_progress_checkTrue by default, the Newton progress check allows the Newton loop to exit if half the allowed iterations have elapsed without a reduction in the *initial* residual. In our experience this usually means the Newton steps are going to fail eventually and we could save some time by quitting early.
NewtonSolver* ContinuationSystem::newton_solver [private]A pointer to the underlying Newton solver used by the DiffSystem. From this pointer, we can get access to all the parameters and options which are available to the 'normal' Newton solver.
unsigned int ContinuationSystem::newton_step [private]Loop counter for nonlinear (Newton) iteration loop.
Real ContinuationSystem::newton_stepgrowth_aggressivenessA measure of how rapidly one should attempt to grow the arclength stepsize based on the number of Newton iterations required to solve the problem. Default value is 1.0, if set to zero, will not try to grow or shrink the arclength stepsize based on the number of Newton iterations required.
Real FEMSystem::numerical_jacobian_h [inherited]If calculating numeric jacobians is required, the FEMSystem will perturb each solution vector entry by numerical_jacobian_h when calculating finite differences.
Real ContinuationSystem::old_continuation_parameterThe system also keeps track of the old value of the continuation parameter.
bool DifferentiableSystem::postprocess_sides [inherited]If postprocess_sides is true (it is false by default), the postprocessing loop will loop over all sides as well as all elements.
Predictor ContinuationSystem::predictor
Real ContinuationSystem::previous_dlambda_ds [private]The old parameter tangent value.
Real ContinuationSystem::previous_ds [private]The previous arcstep length used.
NumericVector<Number>* ContinuationSystem::previous_du_ds [private]The value of du_ds from the previous solution
NumericVector<Number>* ContinuationSystem::previous_u [private]The solution at the previous value of the continuation variable.
bool DifferentiableSystem::print_element_jacobians [inherited]Set print_element_jacobians to true to print each J_elem contribution.
bool DifferentiableSystem::print_jacobian_norms [inherited]Set print_jacobian_norms to true to print |J| whenever it is assembled.
bool DifferentiableSystem::print_jacobians [inherited]Set print_jacobians to true to print J whenever it is assembled.
bool DifferentiableSystem::print_residual_norms [inherited]Set print_residual_norms to true to print |F| whenever it is assembled.
bool DifferentiableSystem::print_residuals [inherited]Set print_residuals to true to print F whenever it is assembled.
bool DifferentiableSystem::print_solution_norms [inherited]Set print_residual_norms to true to print |U| whenever it is used in an assembly() call
bool DifferentiableSystem::print_solutions [inherited]Set print_solutions to true to print U whenever it is used in an assembly() call
Number System::qoi [inherited]Value of the quantity of interest
bool ContinuationSystem::quietIf quiet==false, the System prints extra information about what it is doing.
NumericVector<Number>* ExplicitSystem::rhs [inherited]The system matrix. Implicit systems are characterized by the need to solve the linear system Ax=b. This is the right-hand-side vector b.
RHS_Mode ContinuationSystem::rhs_mode [protected]
AutoPtr<NumericVector<Number> > System::solution [inherited]Data structure to hold solution values.
Real ContinuationSystem::solution_toleranceHow tightly should the Newton iterations attempt to converge ||delta_u|| Defaults to 1.e-6.
bool ContinuationSystem::tangent_initialized [private]False until initialize_tangent() is called
Real ContinuationSystem::ThetaArclength normalization parameter. Defaults to 1.0 (no normalization). Used to ensure that one term in the arclength contstraint equation does not wash out all the others.
Real ContinuationSystem::Theta_LOCAAnother normalization parameter, which is described in the LOCA manual. This one is designed to maintain a relatively 'fixed' value of d(lambda)/ds. It is initially 1.0 and is updated after each solve.
Real DifferentiableSystem::time [inherited]For time-dependent problems, this is the time t at the beginning of the current timestep. But do *not* access this time during an assembly! Use the DiffContext::time value instead to get correct results.
AutoPtr<TimeSolver> DifferentiableSystem::time_solver [inherited]A pointer to the solver object we're going to use. This must be instantiated by the user before solving!
bool DifferentiableSystem::use_fixed_solution [inherited]A boolean to be set to true by systems using elem_fixed_solution, so that the library code will create it. False by default. Warning: if this variable is set to true, it must be before init_data() is called.
Real FEMSystem::verify_analytic_jacobians [inherited]If verify_analytic_jacobian is equal to zero (as it is by default), no numeric jacobians will be calculated unless an overloaded element_time_derivative(), element_constraint(), side_time_derivative(), or side_constraint() function cannot provide an analytic jacobian upon request.
NumericVector<Number>* ContinuationSystem::y [private]Temporary vector 'y' ... stores -du/dlambda, the solution of Ay=G_{}.
NumericVector<Number>* ContinuationSystem::y_old [private]Temporary vector 'y_old' ... stores the previous value of -du/dlambda, which is the solution of Ay=G_{}.
NumericVector<Number>* ContinuationSystem::z [private]Temporary vector 'z' ... the solution of Az = -G
Author

This document was created by man2html, using the manual pages.
Time: 21:43:45 GMT, April 16, 2011