NewmarkSystem Class Reference

#include <newmark_system.h>

Inheritance diagram for NewmarkSystem:

[legend]

List of all members.

Public Types
typedef NewmarkSystem	sys_type
typedef ImplicitSystem	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
	NewmarkSystem (EquationSystems &es, const std::string &name, const unsigned int number)
	~NewmarkSystem ()
virtual void	clear ()
virtual void	reinit ()
virtual void	assemble ()
virtual std::string	system_type () const
void	initial_conditions ()
void	compute_matrix ()
void	update_rhs ()
void	update_u_v_a ()
void	set_newmark_parameters (const Real delta_T=_default_timestep, const Real alpha=_default_alpha, const Real delta=_default_delta)
sys_type &	system ()
virtual void	solve ()
virtual LinearSolver< Number > *	get_linear_solver () const
virtual void	release_linear_solver (LinearSolver< Number > *) const
virtual void	assembly (bool get_residual, bool get_jacobian)
unsigned int	n_linear_iterations () const
Real	final_linear_residual () const
void	attach_shell_matrix (ShellMatrix< Number > *shell_matrix)
void	detach_shell_matrix (void)
ShellMatrix< Number > *	get_shell_matrix (void)
virtual std::pair< unsigned int, Real >	get_linear_solve_parameters () const
virtual void	assemble_residual_derivatives (const ParameterVector &parameters)
virtual std::pair< unsigned int, Real >	sensitivity_solve (const ParameterVector &parameters)
virtual std::pair< unsigned int, Real >	weighted_sensitivity_solve (const ParameterVector &parameters, const ParameterVector &weights)
virtual std::pair< unsigned int, Real >	adjoint_solve (const QoISet &qoi_indices=QoISet())
virtual std::pair< unsigned int, Real >	weighted_sensitivity_adjoint_solve (const ParameterVector &parameters, const ParameterVector &weights, const QoISet &qoi_indices=QoISet())
virtual void	adjoint_qoi_parameter_sensitivity (const QoISet &qoi_indices, const ParameterVector &parameters, SensitivityData &sensitivities)
virtual void	forward_qoi_parameter_sensitivity (const QoISet &qoi_indices, const ParameterVector &parameters, SensitivityData &sensitivities)
virtual void	qoi_parameter_hessian_vector_product (const QoISet &qoi_indices, const ParameterVector &parameters, const ParameterVector &vector, SensitivityData &product)
SparseMatrix< Number > &	add_matrix (const std::string &mat_name)
bool	have_matrix (const std::string &mat_name) const
const SparseMatrix< Number > *	request_matrix (const std::string &mat_name) const
SparseMatrix< Number > *	request_matrix (const std::string &mat_name)
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
virtual void	assemble_qoi (const QoISet &qoi_indices=QoISet())
virtual void	assemble_qoi_derivative (const QoISet &qoi_indices=QoISet())
void	init ()
virtual void	update ()
virtual void	qoi_parameter_sensitivity (const QoISet &qoi_indices, const ParameterVector &parameters, SensitivityData &sensitivities)
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 > *	request_vector (const std::string &vec_name) const
NumericVector< Number > *	request_vector (const std::string &vec_name)
const NumericVector< Number > *	request_vector (const unsigned int vec_num) const
NumericVector< Number > *	request_vector (const unsigned int vec_num)
const NumericVector< Number > &	get_vector (const std::string &vec_name) const
NumericVector< Number > &	get_vector (const std::string &vec_name)
const NumericVector< Number > &	get_vector (const unsigned int vec_num) const
NumericVector< Number > &	get_vector (const unsigned int vec_num)
const std::string &	vector_name (const unsigned int vec_num)
NumericVector< Number > &	add_adjoint_solution (unsigned int i=0)
NumericVector< Number > &	get_adjoint_solution (unsigned int i=0)
const NumericVector< Number > &	get_adjoint_solution (unsigned int i=0) const
NumericVector< Number > &	add_sensitivity_solution (unsigned int i=0)
NumericVector< Number > &	get_sensitivity_solution (unsigned int i=0)
const NumericVector< Number > &	get_sensitivity_solution (unsigned int i=0) const
NumericVector< Number > &	add_weighted_sensitivity_adjoint_solution (unsigned int i=0)
NumericVector< Number > &	get_weighted_sensitivity_adjoint_solution (unsigned int i=0)
const NumericVector< Number > &	get_weighted_sensitivity_adjoint_solution (unsigned int i=0) const
NumericVector< Number > &	add_weighted_sensitivity_solution ()
NumericVector< Number > &	get_weighted_sensitivity_solution ()
const NumericVector< Number > &	get_weighted_sensitivity_solution () const
NumericVector< Number > &	add_adjoint_rhs (unsigned int i=0)
NumericVector< Number > &	get_adjoint_rhs (unsigned int i=0)
const NumericVector< Number > &	get_adjoint_rhs (unsigned int i=0) const
NumericVector< Number > &	add_sensitivity_rhs (unsigned int i=0)
NumericVector< Number > &	get_sensitivity_rhs (unsigned int i=0)
const NumericVector< Number > &	get_sensitivity_rhs (unsigned int i=0) const
unsigned int	n_vectors () const
unsigned int	n_vars () const
unsigned int	n_dofs () const
unsigned int	n_active_dofs () const
unsigned int	n_constrained_dofs () const
unsigned int	n_local_dofs () const
unsigned int	add_variable (const std::string &var, const FEType &type, const std::set< subdomain_id_type > *const active_subdomains=NULL)
unsigned int	add_variable (const std::string &var, const Order order=FIRST, const FEFamily=LAGRANGE, const std::set< subdomain_id_type > *const active_subdomains=NULL)
const Variable &	variable (unsigned int var) const
bool	has_variable (const std::string &var) const
const std::string &	variable_name (const unsigned int i) const
unsigned short int	variable_number (const std::string &var) const
const FEType &	variable_type (const unsigned int i) const
const FEType &	variable_type (const std::string &var) const
Real	calculate_norm (NumericVector< Number > &v, unsigned int var=0, FEMNormType norm_type=L2) const
Real	calculate_norm (NumericVector< Number > &v, const SystemNorm &norm) const
void	read_header (Xdr &io, const std::string &version, const bool read_header=true, const bool read_additional_data=true, const bool read_legacy_format=false)
void	read_legacy_data (Xdr &io, const bool read_additional_data=true)
void	read_serialized_data (Xdr &io, const bool read_additional_data=true)
void	read_parallel_data (Xdr &io, const bool read_additional_data)
void	write_header (Xdr &io, const std::string &version, const bool write_additional_data) const
void	write_serialized_data (Xdr &io, const bool write_additional_data=true) const
void	write_parallel_data (Xdr &io, const bool write_additional_data) const
std::string	get_info () const
void	attach_init_function (void fptr(EquationSystems &es, const std::string &name))
void	attach_assemble_function (void fptr(EquationSystems &es, const std::string &name))
void	attach_constraint_function (void fptr(EquationSystems &es, const std::string &name))
void	attach_QOI_function (void fptr(EquationSystems &es, const std::string &name, const QoISet &qoi_indices))
void	attach_QOI_derivative (void fptr(EquationSystems &es, const std::string &name, const QoISet &qoi_indices))
virtual void	user_initialization ()
virtual void	user_assembly ()
virtual void	user_constrain ()
virtual void	user_QOI (const QoISet &qoi_indices)
virtual void	user_QOI_derivative (const QoISet &qoi_indices)
virtual void	re_update ()
virtual void	restrict_vectors ()
virtual void	prolong_vectors ()
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
AutoPtr< LinearSolver< Number > >	linear_solver
SparseMatrix< Number > *	matrix
NumericVector< Number > *	rhs
bool	assemble_before_solve
AutoPtr< NumericVector< Number > >	solution
AutoPtr< NumericVector< Number > >	current_local_solution
std::vector< Number >	qoi
Protected Types
typedef std::map< std::string, std::pair< unsigned int, unsigned int > >	Counts
Protected Member Functions
virtual void	init_data ()
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
unsigned int	_n_linear_iterations
Real	_final_linear_residual
ShellMatrix< Number > *	_shell_matrix
Static Protected Attributes
static Counts	_counts
static Threads::atomic < unsigned int >	_n_objects
static Threads::spin_mutex	_mutex
Private Attributes
Real	_a_0
Real	_a_1
Real	_a_2
Real	_a_3
Real	_a_4
Real	_a_5
Real	_a_6
Real	_a_7
bool	_finished_assemble
Static Private Attributes
static const Real	_default_alpha = .25
static const Real	_default_delta = .5
static const Real	_default_timestep = 1.

---

Detailed Description

This class contains a specific system class. It provides an implicit time integration scheme known as the Newmark method.

In the algorithm implemented here the system is solved for displacements. Curently the Newmark scheme is implemented for constant time step sizes only. This time step is stored in the EquationSystems parameter named "Newmark \p time \p step". For the case of constant time steps the matrix only has to be assembled once, whereas the rhs has to be updated in each timestep. Default values of the Newmark parameters alpha and delta used for time integration are provided. For details refer to the examples section.

Definition at line 53 of file newmark_system.h.

---

Member Typedef Documentation

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

Definition at line 256 of file implicit_system.h.

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

Definition at line 405 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 255 of file implicit_system.h.

typedef ImplicitSystem LinearImplicitSystem::Parent [inherited]

The type of the parent.

Reimplemented from ImplicitSystem.

Definition at line 72 of file linear_implicit_system.h.

typedef NewmarkSystem NewmarkSystem::sys_type

The type of system.

Reimplemented from LinearImplicitSystem.

Definition at line 74 of file newmark_system.h.

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

Vector iterator typedefs.

Definition at line 404 of file system.h.

---

Constructor & Destructor Documentation

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

Constructor. Optionally initializes required data structures.

Definition at line 44 of file newmark_system.C.

References _default_alpha, _default_delta, _default_timestep, ImplicitSystem::add_matrix(), System::add_vector(), EquationSystems::parameters, and Parameters::set().

                                                         :  
  LinearImplicitSystem (es, name, number),
  _a_0                 (1./(_default_alpha*_default_timestep*_default_timestep)),
  _a_1                 (_default_delta/(_default_alpha*_default_timestep)),
  _a_2                 (1./(_default_alpha*_default_timestep)),
  _a_3                 (1./(2.*_default_alpha)-1.),
  _a_4                 (_default_delta/_default_alpha -1.),
  _a_5                 (_default_timestep/2.*(_default_delta/_default_alpha-2.)),
  _a_6                 (_default_timestep*(1.-_default_delta)),
  _a_7                 (_default_delta*_default_timestep),
  _finished_assemble   (false)
  
{
  // default values of the newmark parameters
  es.parameters.set<Real>("Newmark alpha") = _default_alpha;
  es.parameters.set<Real>("Newmark delta") = _default_delta;

  // time step size.
  // should be handled at a later stage through EquationSystems?
  es.parameters.set<Real>("Newmark time step") = _default_timestep;

  // add additional matrices and vectors that will be used in the
  // newmark algorithm to the data structure
  // functions LinearImplicitSystem::add_matrix and LinearImplicitSystem::add_vector
  // are used so we do not have to bother about initialization and
  // dof mapping

  // system matrices
  this->add_matrix ("stiffness");
  this->add_matrix ("damping");
  this->add_matrix ("mass");

  // load vector  
  this->add_vector ("force");

  // the displacement and the time derivatives  
  this->add_vector ("displacement");
  this->add_vector ("velocity");
  this->add_vector ("acceleration");

  // contributions to the rhs
  this->add_vector ("rhs_m");
  this->add_vector ("rhs_c");
  
  // results from the previous time step
  this->add_vector ("old_solution");
  this->add_vector ("old_acceleration");
}

NewmarkSystem::~NewmarkSystem

(

)

Destructor.

Definition at line 97 of file newmark_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 1374 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 1366 of file system.h.

References System::_active.

{
  return _active;  
}

NumericVector< Number > & System::add_adjoint_rhs

(

unsigned int

i = 0

)

[inherited]

Returns:

a reference to one of the system's adjoint rhs vectors, by default the one corresponding to the first qoi. Creates the vector if it doesn't already exist.

Definition at line 824 of file system.C.

References System::add_vector().

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

{
  OStringStream adjoint_rhs_name;
  adjoint_rhs_name << "adjoint_rhs" << i;

  return this->add_vector(adjoint_rhs_name.str());
}

NumericVector< Number > & System::add_adjoint_solution

(

unsigned int

i = 0

)

[inherited]

Returns:

a reference to one of the system's adjoint solution vectors, by default the one corresponding to the first qoi. Creates the vector if it doesn't already exist.

Definition at line 764 of file system.C.

References System::add_vector().

Referenced by ImplicitSystem::adjoint_solve().

{
  OStringStream adjoint_name;
  adjoint_name << "adjoint_solution" << i;

  return this->add_vector(adjoint_name.str());
}

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 175 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().

{
  // 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;
}

NumericVector< Number > & System::add_sensitivity_rhs

(

unsigned int

i = 0

)

[inherited]

Returns:

a reference to one of the system's sensitivity rhs vectors, by default the one corresponding to the first parameter. Creates the vector if it doesn't already exist.

Definition at line 854 of file system.C.

References System::add_vector().

Referenced by ImplicitSystem::assemble_residual_derivatives().

{
  OStringStream sensitivity_rhs_name;
  sensitivity_rhs_name << "sensitivity_rhs" << i;

  return this->add_vector(sensitivity_rhs_name.str());
}

NumericVector< Number > & System::add_sensitivity_solution

(

unsigned int

i = 0

)

[inherited]

Returns:

a reference to one of the system's solution sensitivity vectors, by default the one corresponding to the first parameter. Creates the vector if it doesn't already exist.

Definition at line 714 of file system.C.

References System::add_vector().

{
  OStringStream sensitivity_name;
  sensitivity_name << "sensitivity_solution" << i;

  return this->add_vector(sensitivity_name.str());
}

unsigned int System::add_variable	(	const std::string &	var,
		const Order	order = FIRST,
		const FEFamily	family = LAGRANGE,
		const std::set< subdomain_id_type > *const	active_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 923 of file system.C.

References System::add_variable().

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

unsigned int System::add_variable	(	const std::string &	var,
		const FEType &	type,
		const std::set< subdomain_id_type > *const	active_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 884 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;
}

NumericVector< Number > & System::add_vector	(	const std::string &	vec_name,
		const bool	projections = 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 549 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 System::add_adjoint_rhs(), System::add_adjoint_solution(), System::add_sensitivity_rhs(), System::add_sensitivity_solution(), ExplicitSystem::add_system_rhs(), System::add_weighted_sensitivity_solution(), UnsteadySolver::init(), ContinuationSystem::init_data(), 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;
}

NumericVector< Number > & System::add_weighted_sensitivity_adjoint_solution

(

unsigned int

i = 0

)

[inherited]

Returns:

a reference to one of the system's weighted sensitivity adjoint solution vectors, by default the one corresponding to the first qoi. Creates the vector if it doesn't already exist.

Definition at line 794 of file system.C.

References System::get_vector().

Referenced by ImplicitSystem::weighted_sensitivity_adjoint_solve().

{
  OStringStream adjoint_name;
  adjoint_name << "weighted_sensitivity_adjoint_solution" << i;

  return this->get_vector(adjoint_name.str());
}

NumericVector< Number > & System::add_weighted_sensitivity_solution

(

)

[inherited]

Returns:

a reference to the solution of the last weighted sensitivity solve Creates the vector if it doesn't already exist.

Definition at line 743 of file system.C.

References System::add_vector().

Referenced by ImplicitSystem::weighted_sensitivity_solve().

{
  return this->add_vector("weighted_sensitivity_solution");
}

void ImplicitSystem::adjoint_qoi_parameter_sensitivity	(	const QoISet &	qoi_indices,
		const ParameterVector &	parameters,
		SensitivityData &	sensitivities
	)			[virtual, inherited]

Solves for the derivative of each of the system's quantities of interest q in qoi[qoi_indices] with respect to each parameter in parameters, placing the result for qoi i and parameter j into sensitivities[i][j].

Uses adjoint_solve() and the adjoint sensitivity method.

Currently uses finite differenced derivatives (partial q / partial p) and (partial R / partial p).

Reimplemented from System.

Definition at line 657 of file implicit_system.C.

References SensitivityData::allocate_data(), QoISet::has_index(), and ParameterVector::size().

{
  const unsigned int Np = parameters.size();
  const unsigned int Nq = qoi.size();

  // We currently get partial derivatives via central differencing
  const Real delta_p = TOLERANCE;

  // 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(qoi_indices);

  // Get ready to fill in senstivities:
  sensitivities.allocate_data(qoi_indices, *this, parameters);

  // 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 j=0; j != Np; ++j)
    {
      // (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[j];
      // Number old_qoi = this->qoi;

      *parameters[j] = old_parameter - delta_p;
      this->assemble_qoi(qoi_indices);
      std::vector<Number> qoi_minus = this->qoi;

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

      *parameters[j] = old_parameter + delta_p;
      this->assemble_qoi(qoi_indices);
      std::vector<Number>& qoi_plus = this->qoi;

      std::vector<Number> partialq_partialp(Nq, 0);
      for (unsigned int i=0; i != Nq; ++i)
        if (qoi_indices.has_index(i))
          partialq_partialp[i] = (qoi_plus[i] - qoi_minus[i]) / (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[j] = old_parameter;

      for (unsigned int i=0; i != Nq; ++i)
        if (qoi_indices.has_index(i))
          sensitivities[i][j] = partialq_partialp[i] -
            partialR_partialp->dot(this->get_adjoint_solution(i));
    }

  // All parameters have been reset.
  // We didn't cache the original rhs or matrix for memory reasons,
  // but we can restore them to a state consistent solution -
  // principle of least surprise.
  this->assembly(true, true);
  this->rhs->close();
  this->matrix->close();
  this->assemble_qoi(qoi_indices);
}

std::pair< unsigned int, Real > ImplicitSystem::adjoint_solve

(

const QoISet &

qoi_indices = QoISet()

)

[virtual, inherited]

Assembles & solves the linear system (dR/du)^T*z = dq/du, for those quantities of interest q specified by qoi_indices.

Leave qoi_indices empty to solve all adjoint problems.

Returns a pair with the total number of linear iterations performed and the (sum of the) final residual norms

Reimplemented from System.

Definition at line 342 of file implicit_system.C.

References System::add_adjoint_solution(), System::assemble_before_solve, ExplicitSystem::assemble_qoi_derivative(), ImplicitSystem::assembly(), DofMap::enforce_constraints_exactly(), System::get_adjoint_rhs(), System::get_adjoint_solution(), System::get_dof_map(), ImplicitSystem::get_linear_solve_parameters(), ImplicitSystem::get_linear_solver(), QoISet::has_index(), ImplicitSystem::matrix, System::qoi, ImplicitSystem::release_linear_solver(), and LinearSolver< T >::solve().

{
  // Log how long the linear solve takes.
  START_LOG("adjoint_solve()", "ImplicitSystem");

  // The forward system should now already be solved.
  // Now assemble it's adjoint.

  if (this->assemble_before_solve)
    {
      // Build the Jacobian
      this->assembly(false, true);
      this->matrix->close();

      // Take the discrete adjoint
      matrix->get_transpose(*matrix);

      // Reset and build the RHS from the QOI derivative
      this->assemble_qoi_derivative(qoi_indices);
    }

  // The adjoint problem is linear
  LinearSolver<Number> *linear_solver = this->get_linear_solver();

  // Our iteration counts and residuals will be sums of the individual
  // results
  std::pair<unsigned int, Real> solver_params =
    this->get_linear_solve_parameters();
  std::pair<unsigned int, Real> totalrval = std::make_pair(0,0.0);

  for (unsigned int i=0; i != this->qoi.size(); ++i)
    if (qoi_indices.has_index(i))
      {
        const std::pair<unsigned int, Real> rval =
          linear_solver->solve (*matrix, this->add_adjoint_solution(i),
                                 this->get_adjoint_rhs(i),
                                 solver_params.second,
                                 solver_params.first);

        totalrval.first  += rval.first;
        totalrval.second += rval.second;
      }

  this->release_linear_solver(linear_solver);

  // The linear solver may not have fit our constraints exactly
#ifdef LIBMESH_ENABLE_AMR
  for (unsigned int i=0; i != this->qoi.size(); ++i)
    if (qoi_indices.has_index(i))
      this->get_dof_map().enforce_constraints_exactly
        (*this, &this->get_adjoint_solution(i));
#endif

  // Stop logging the nonlinear solve
  STOP_LOG("adjoint_solve()", "ImplicitSystem");

  return totalrval;
}

void NewmarkSystem::assemble

(

)

[virtual]

Assemble the linear system. Does not actually call the solver.

Reimplemented from LinearImplicitSystem.

Definition at line 137 of file newmark_system.C.

References _finished_assemble, compute_matrix(), and initial_conditions().

{
  if (!_finished_assemble)
    {
      // prepare matrix with the help of the _dof_map, 
      // fill with sparsity pattern
      LinearImplicitSystem::assemble();

      // compute the effective system matrix
      this->compute_matrix();

      // apply initial conditions
      this->initial_conditions();
      
      _finished_assemble = true;
    }
}

void ExplicitSystem::assemble_qoi

(

const QoISet &

qoi_indices = QoISet()

)

[virtual, inherited]

Prepares qoi for quantity of interest assembly, then calls user qoi function. Can be overloaded in derived classes.

Reimplemented from System.

Reimplemented in FEMSystem.

Definition at line 89 of file explicit_system.C.

References System::assemble_qoi(), QoISet::has_index(), and System::qoi.

{
  // The user quantity of interest assembly gets to expect to
  // accumulate on initially zero values
  for (unsigned int i=0; i != qoi.size(); ++i)
    if (qoi_indices.has_index(i))
      qoi[i] = 0;

  Parent::assemble_qoi (qoi_indices);
}

void ExplicitSystem::assemble_qoi_derivative

(

const QoISet &

qoi_indices = QoISet()

)

[virtual, inherited]

Prepares adjoint_rhs for quantity of interest derivative assembly, then calls user qoi derivative function. Can be overloaded in derived classes.

Reimplemented from System.

Reimplemented in FEMSystem.

Definition at line 102 of file explicit_system.C.

References System::add_adjoint_rhs(), System::assemble_qoi_derivative(), QoISet::has_index(), and System::qoi.

Referenced by ImplicitSystem::adjoint_solve(), and ImplicitSystem::weighted_sensitivity_adjoint_solve().

{
  // The user quantity of interest derivative assembly gets to expect
  // to accumulate on initially zero vectors
  for (unsigned int i=0; i != qoi.size(); ++i)
    if (qoi_indices.has_index(i))
      this->add_adjoint_rhs(i).zero();

  Parent::assemble_qoi_derivative (qoi_indices);
}

void ImplicitSystem::assemble_residual_derivatives

(

const ParameterVector &

parameters

)

[virtual, inherited]

Residual parameter derivative function.

Uses finite differences by default.

This will assemble the sensitivity rhs vectors to hold -(partial R / partial p_i), making them ready to solve the forward sensitivity equation.

Can be overloaded in derived classes.

Reimplemented from System.

Definition at line 622 of file implicit_system.C.

References System::add_sensitivity_rhs(), ImplicitSystem::assembly(), NumericVector< T >::close(), ExplicitSystem::rhs, and ParameterVector::size().

Referenced by ImplicitSystem::sensitivity_solve().

{
  const unsigned int Np = parameters.size();
  Real deltap = TOLERANCE;

  for (unsigned int p=0; p != Np; ++p)
    {
      NumericVector<Number> &sensitivity_rhs = this->add_sensitivity_rhs(p);

      // Approximate -(partial R / partial p) by
      // (R(p-dp) - R(p+dp)) / (2*dp)

      Number old_parameter = *parameters[p];
      *parameters[p] -= deltap;

      this->assembly(true, false);
      this->rhs->close();
      sensitivity_rhs = *this->rhs;

      *parameters[p] = old_parameter + deltap;

      this->assembly(true, false);
      this->rhs->close();

      sensitivity_rhs -= *this->rhs;
      sensitivity_rhs /= (2*deltap);
      sensitivity_rhs.close();

      *parameters[p] = old_parameter;
    }
}

void LinearImplicitSystem::assembly	(	bool	get_residual,
		bool	get_jacobian
	)			[virtual, inherited]

Assembles a residual in rhs and/or a jacobian in matrix, as requested.

Reimplemented from ImplicitSystem.

Definition at line 333 of file linear_implicit_system.C.

References LinearImplicitSystem::assemble(), ImplicitSystem::matrix, ExplicitSystem::rhs, and System::solution.

{
  // Residual R(u(p),p) := A(p)*u(p) - b(p)
  // partial R / partial u = A

  this->assemble();
  this->rhs->close();
  this->matrix->close();

  *(this->rhs) *= -1.0;
  this->rhs->add_vector(*this->solution, *this->matrix);
}

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 1316 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 1326 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 1306 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, const QoISet &qoi_indices

)

[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, const QoISet &qoi_indices

)

[inherited]

Register a user function for evaluating the quantities of interest, whose values should be placed in System::qoi

void LinearImplicitSystem::attach_shell_matrix

(

ShellMatrix< Number > *

shell_matrix

)

[inherited]

This function enables the user to provide a shell matrix, i.e. a matrix that is not stored element-wise, but as a function. When you register your shell matrix using this function, calling solve() will no longer use the matrix member but the registered shell matrix instead. You can reset this behaviour to its original state by supplying a NULL pointer to this function.

Definition at line 125 of file linear_implicit_system.C.

References LinearImplicitSystem::_shell_matrix.

Referenced by LinearImplicitSystem::detach_shell_matrix().

{
  _shell_matrix = shell_matrix;
}

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

Real System::calculate_norm	(	NumericVector< Number > &	v,
		unsigned int	var = 0,
		FEMNormType	norm_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 1056 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;
}

void NewmarkSystem::clear

(

)

[virtual]

Clear all the data structures associated with the system.

Reimplemented from LinearImplicitSystem.

Definition at line 104 of file newmark_system.C.

References _default_alpha, _default_delta, _default_timestep, _finished_assemble, System::get_equation_systems(), EquationSystems::parameters, and Parameters::set().

Referenced by ~NewmarkSystem().

{
  // use parent clear this will also clear the
  // matrices and vectors added in the constructor
  LinearImplicitSystem::clear();

  // Get a reference to the EquationSystems
  EquationSystems& es =
    this->get_equation_systems();
  
  // default values of the newmark parameters
  es.parameters.set<Real>("Newmark alpha") = _default_alpha;
  es.parameters.set<Real>("Newmark delta") = _default_delta;

  // time step size.  should be handled at a later stage through EquationSystems?
  es.parameters.set<Real>("Newmark time step") = _default_timestep;

  // set bool to false
  _finished_assemble = false;
}

bool System::compare	(	const System &	other_system,
		const Real	threshold,
		const bool	verbose
	)			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 399 of file system.C.

References System::_can_add_vectors, System::_sys_name, System::_vectors, 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 NewmarkSystem::compute_matrix

(

)

Compute the global matrix by adding up scaled mass damping and stiffness matrix.

Definition at line 176 of file newmark_system.C.

References _a_0, _a_1, ImplicitSystem::get_matrix(), and ImplicitSystem::matrix.

Referenced by assemble().

{
  // close the component matrices
  this->get_matrix ("stiffness").close();
  this->get_matrix ("mass"     ).close();
  this->get_matrix ("damping"  ).close();

  // close & zero the system matrix
  this->matrix->close (); this->matrix->zero();

  // add up the matrices
  this->matrix->add (1.,   this->get_matrix ("stiffness"));
  this->matrix->add (_a_0, this->get_matrix ("mass"));
  this->matrix->add (_a_1, this->get_matrix ("damping"));

}

Number System::current_solution

(

const unsigned int

global_dof_number

)

const [inherited]

Returns:

the current solution for the specified global DOF.

Definition at line 117 of file system.C.

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

Referenced by ExactSolution::_compute_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 1382 of file system.h.

References System::_active.

{
  _active = false;
}

void LinearImplicitSystem::detach_shell_matrix

(

void

)

[inline, inherited]

Detaches a shell matrix. Same as attach_shell_matrix(NULL).

Definition at line 161 of file linear_implicit_system.h.

References LinearImplicitSystem::attach_shell_matrix().

{ attach_shell_matrix(NULL); }

Real LinearImplicitSystem::final_linear_residual

(

)

const [inline, inherited]

Returns the final residual for the linear system solve.

Definition at line 145 of file linear_implicit_system.h.

References LinearImplicitSystem::_final_linear_residual.

{ return _final_linear_residual; }

void ImplicitSystem::forward_qoi_parameter_sensitivity	(	const QoISet &	qoi_indices,
		const ParameterVector &	parameters,
		SensitivityData &	sensitivities
	)			[virtual, inherited]

Solves for the derivative of each of the system's quantities of interest q in qoi[qoi_indices] with respect to each parameter in parameters, placing the result for qoi i and parameter j into sensitivities[i][j].

Uses the forward sensitivity method.

Currently uses finite differenced derivatives (partial q / partial p) and (partial R / partial p).

Reimplemented from System.

Definition at line 749 of file implicit_system.C.

References SensitivityData::allocate_data(), QoISet::has_index(), and ParameterVector::size().

{
  const unsigned int Np = parameters.size();
  const unsigned int Nq = qoi.size();

  // We currently get partial derivatives via central differencing
  const Real delta_p = TOLERANCE;

  // 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 solve for (partial u / partial p) for each parameter:
  // J * (partial u / partial p) = - (partial R / partial p)

  this->sensitivity_solve(parameters);

  // Get ready to fill in senstivities:
  sensitivities.allocate_data(qoi_indices, *this, parameters);

  // We use the identity:
  // dq/dp = (partial q / partial p) + (partial q / partial u) *
  //         (partial u / partial p)
 
  // We get (partial q / partial u) from the user
  this->assemble_qoi_derivative(qoi_indices);

  for (unsigned int j=0; j != Np; ++j)
    {
      // (partial q / partial p) ~= (q(p+dp)-q(p-dp))/(2*dp)

      Number old_parameter = *parameters[j];

      *parameters[j] = old_parameter - delta_p;
      this->assemble_qoi();
      std::vector<Number> qoi_minus = this->qoi;

      *parameters[j] = old_parameter + delta_p;
      this->assemble_qoi();
      std::vector<Number>& qoi_plus = this->qoi;

      std::vector<Number> partialq_partialp(Nq, 0);
      for (unsigned int i=0; i != Nq; ++i)
        if (qoi_indices.has_index(i))
          partialq_partialp[i] = (qoi_plus[i] - qoi_minus[i]) / (2.*delta_p);

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

      for (unsigned int i=0; i != Nq; ++i)
        if (qoi_indices.has_index(i))
          sensitivities[i][j] = partialq_partialp[i] +
            this->get_adjoint_rhs(i).dot(this->get_sensitivity_solution(i));
    }

  // All parameters have been reset.
  // We didn't cache the original rhs or matrix for memory reasons,
  // but we can restore them to a state consistent solution -
  // principle of least surprise.
  this->assembly(true, true);
  this->rhs->close();
  this->matrix->close();
  this->assemble_qoi(qoi_indices);
}

const NumericVector< Number > & System::get_adjoint_rhs

(

unsigned int

i = 0

)

const [inherited]

Returns:

a reference to one of the system's adjoint rhs vectors, by default the one corresponding to the first qoi.

Definition at line 844 of file system.C.

References System::get_vector().

{
  OStringStream adjoint_rhs_name;
  adjoint_rhs_name << "adjoint_rhs" << i;

  return this->get_vector(adjoint_rhs_name.str());
}

NumericVector< Number > & System::get_adjoint_rhs

(

unsigned int

i = 0

)

[inherited]

Returns:

a reference to one of the system's adjoint rhs vectors, by default the one corresponding to the first qoi. This what the user's QoI derivative code should assemble when setting up an adjoint problem

Definition at line 834 of file system.C.

References System::get_vector().

Referenced by ImplicitSystem::adjoint_solve(), FEMSystem::assemble_qoi_derivative(), and ImplicitSystem::weighted_sensitivity_adjoint_solve().

{
  OStringStream adjoint_rhs_name;
  adjoint_rhs_name << "adjoint_rhs" << i;

  return this->get_vector(adjoint_rhs_name.str());
}

const NumericVector< Number > & System::get_adjoint_solution

(

unsigned int

i = 0

)

const [inherited]

Returns:

a reference to one of the system's adjoint solution vectors, by default the one corresponding to the first qoi.

Definition at line 784 of file system.C.

References System::get_vector().

{
  OStringStream adjoint_name;
  adjoint_name << "adjoint_solution" << i;

  return this->get_vector(adjoint_name.str());
}

NumericVector< Number > & System::get_adjoint_solution

(

unsigned int

i = 0

)

[inherited]

Returns:

a reference to one of the system's adjoint solution vectors, by default the one corresponding to the first qoi.

Definition at line 774 of file system.C.

References System::get_vector().

Referenced by ImplicitSystem::adjoint_solve(), AdjointResidualErrorEstimator::estimate_error(), and ImplicitSystem::weighted_sensitivity_adjoint_solve().

{
  OStringStream adjoint_name;
  adjoint_name << "adjoint_solution" << i;

  return this->get_vector(adjoint_name.str());
}

DofMap & System::get_dof_map

(

)

[inline, inherited]

Returns:

a writeable reference to this system's _dof_map.

Definition at line 1358 of file system.h.

References System::_dof_map.

{
  return *_dof_map;
}

const DofMap & System::get_dof_map

(

)

const [inline, inherited]

Returns:

a constant reference to this system's _dof_map.

Definition at line 1350 of file system.h.

References System::_dof_map.

Referenced by __libmesh_petsc_diff_solver_jacobian(), __libmesh_petsc_diff_solver_residual(), ExactSolution::_compute_error(), HPCoarsenTest::add_projection(), ImplicitSystem::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(), ImplicitSystem::sensitivity_solve(), UnsteadySolver::solve(), PetscDiffSolver::solve(), NewtonSolver::solve(), ImplicitSystem::weighted_sensitivity_adjoint_solve(), ImplicitSystem::weighted_sensitivity_solve(), EnsightIO::write_scalar_ascii(), and EnsightIO::write_vector_ascii().

{
  return *_dof_map;
}

EquationSystems& System::get_equation_systems

(

)

[inline, inherited]

Returns:

a reference to this system's parent EquationSystems object.

Definition at line 383 of file system.h.

References System::_equation_systems.

{ return _equation_systems; }

const EquationSystems& System::get_equation_systems

(

)

const [inline, inherited]

Returns:

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

Definition at line 378 of file system.h.

References System::_equation_systems.

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

{ return _equation_systems; }

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 QuadratureRules::name().

Referenced by ReferenceCounter::print_info().

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

  std::ostringstream out;
  
  out << '\n'
      << " ---------------------------------------------------------------------------- \n"
      << "| Reference count information                                                |\n"
      << " ---------------------------------------------------------------------------- \n";
  
  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:\n"
          << "|  Creations:    " << creations    << '\n'
          << "|  Destructions: " << destructions << '\n';
    }
  
  out << " ---------------------------------------------------------------------------- \n";

  return out.str();

#else

  return "";
  
#endif
}

std::string System::get_info

(

)

const [inherited]

Returns:

a string containing information about the system.

Definition at line 1233 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 << "\"\n"
      << "    Type \""  << this->system_type() << "\"\n"
      << "    Variables=";
  
  for (unsigned int vn=0; vn<this->n_vars(); vn++)
      out << "\"" << this->variable_name(vn) << "\" ";
     
  out << '\n';

  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 << '\n' << "    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 << '\n';
      
  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 << '\n';
      
  out << "    n_dofs()="             << this->n_dofs()             << '\n';
  out << "    n_local_dofs()="       << this->n_local_dofs()       << '\n';
#ifdef LIBMESH_ENABLE_AMR
  out << "    n_constrained_dofs()=" << this->n_constrained_dofs() << '\n';
#endif

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

std::pair< unsigned int, Real > ImplicitSystem::get_linear_solve_parameters

(

)

const [virtual, inherited]

Returns an integer corresponding to the upper iteration count limit and a Real corresponding to the convergence tolerance to be used in linear adjoint and/or sensitivity solves

Reimplemented in DifferentiableSystem, and NonlinearImplicitSystem.

Definition at line 1059 of file implicit_system.C.

References System::get_equation_systems().

Referenced by ImplicitSystem::adjoint_solve(), ImplicitSystem::sensitivity_solve(), ImplicitSystem::weighted_sensitivity_adjoint_solve(), and ImplicitSystem::weighted_sensitivity_solve().

{
  return std::make_pair(this->get_equation_systems().parameters.get<unsigned int>("linear solver maximum iterations"),
                        this->get_equation_systems().parameters.get<Real>("linear solver tolerance"));
}

LinearSolver< Number > * LinearImplicitSystem::get_linear_solver

(

)

const [virtual, inherited]

Returns a pointer to a linear solver appropriate for use in adjoint and/or sensitivity solves

Reimplemented from ImplicitSystem.

Definition at line 320 of file linear_implicit_system.C.

References LinearImplicitSystem::linear_solver.

{
  return linear_solver.get();
}

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 245 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 226 of file implicit_system.C.

References ImplicitSystem::_matrices.

Referenced by compute_matrix(), EigenTimeSolver::solve(), and 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);
}

MeshBase & System::get_mesh

(

)

[inline, inherited]

Returns:

a reference to this systems's _mesh.

Definition at line 1342 of file system.h.

References System::_mesh.

{
  return _mesh;
}

const MeshBase & System::get_mesh

(

)

const [inline, inherited]

Returns:

a constant reference to this systems's _mesh.

Definition at line 1334 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;
}

const NumericVector< Number > & System::get_sensitivity_rhs

(

unsigned int

i = 0

)

const [inherited]

Returns:

a reference to one of the system's sensitivity rhs vectors, by default the one corresponding to the first parameter.

Definition at line 874 of file system.C.

References System::get_vector().

{
  OStringStream sensitivity_rhs_name;
  sensitivity_rhs_name << "sensitivity_rhs" << i;

  return this->get_vector(sensitivity_rhs_name.str());
}

NumericVector< Number > & System::get_sensitivity_rhs

(

unsigned int

i = 0

)

[inherited]

Returns:

a reference to one of the system's sensitivity rhs vectors, by default the one corresponding to the first parameter. By default these vectors are built by the library, using finite differences, when assemble_residual_derivatives() is called.

When assembled, this vector should hold -(partial R / partial p_i)

Definition at line 864 of file system.C.

References System::get_vector().

Referenced by ImplicitSystem::sensitivity_solve().

{
  OStringStream sensitivity_rhs_name;
  sensitivity_rhs_name << "sensitivity_rhs" << i;

  return this->get_vector(sensitivity_rhs_name.str());
}

const NumericVector< Number > & System::get_sensitivity_solution

(

unsigned int

i = 0

)

const [inherited]

Returns:

a reference to one of the system's solution sensitivity vectors, by default the one corresponding to the first parameter.

Definition at line 733 of file system.C.

References System::get_vector().

{
  OStringStream sensitivity_name;
  sensitivity_name << "sensitivity_solution" << i;

  return this->get_vector(sensitivity_name.str());
}

NumericVector< Number > & System::get_sensitivity_solution

(

unsigned int

i = 0

)

[inherited]

Returns:

a reference to one of the system's solution sensitivity vectors, by default the one corresponding to the first parameter.

Definition at line 723 of file system.C.

References System::get_vector().

Referenced by ImplicitSystem::sensitivity_solve().

{
  OStringStream sensitivity_name;
  sensitivity_name << "sensitivity_solution" << i;

  return this->get_vector(sensitivity_name.str());
}

ShellMatrix<Number>* LinearImplicitSystem::get_shell_matrix

(

void

)

[inline, inherited]

Returns a pointer to the currently attached shell matrix, if any, or NULL else.

Definition at line 167 of file linear_implicit_system.h.

References LinearImplicitSystem::_shell_matrix.

{ return _shell_matrix; }

NumericVector< Number > & System::get_vector

(

const unsigned int

vec_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 682 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);
}

const NumericVector< Number > & System::get_vector

(

const unsigned int

vec_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 666 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 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 647 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 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 628 of file system.C.

References System::_vectors.

Referenced by System::add_weighted_sensitivity_adjoint_solution(), UnsteadySolver::advance_timestep(), AdaptiveTimeSolver::advance_timestep(), System::compare(), UnsteadySolver::du(), System::get_adjoint_rhs(), System::get_adjoint_solution(), System::get_sensitivity_rhs(), System::get_sensitivity_solution(), System::get_weighted_sensitivity_adjoint_solution(), System::get_weighted_sensitivity_solution(), initial_conditions(), UnsteadySolver::solve(), TwostepTimeSolver::solve(), FrequencySystem::solve(), update_rhs(), and 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);
}

const NumericVector< Number > & System::get_weighted_sensitivity_adjoint_solution

(

unsigned int

i = 0

)

const [inherited]

Returns:

a reference to one of the system's weighted sensitivity adjoint solution vectors, by default the one corresponding to the first qoi.

Definition at line 814 of file system.C.

References System::get_vector().

{
  OStringStream adjoint_name;
  adjoint_name << "weighted_sensitivity_adjoint_solution" << i;

  return this->get_vector(adjoint_name.str());
}

NumericVector< Number > & System::get_weighted_sensitivity_adjoint_solution

(

unsigned int

i = 0

)

[inherited]

Returns:

a reference to one of the system's weighted sensitivity adjoint solution vectors, by default the one corresponding to the first qoi.

Definition at line 804 of file system.C.

References System::get_vector().

Referenced by ImplicitSystem::weighted_sensitivity_adjoint_solve().

{
  OStringStream adjoint_name;
  adjoint_name << "weighted_sensitivity_adjoint_solution" << i;

  return this->get_vector(adjoint_name.str());
}

const NumericVector< Number > & System::get_weighted_sensitivity_solution

(

)

const [inherited]

Returns:

a reference to the solution of the last weighted sensitivity solve

Definition at line 757 of file system.C.

References System::get_vector().

{
  return this->get_vector("weighted_sensitivity_solution");
}

NumericVector< Number > & System::get_weighted_sensitivity_solution

(

)

[inherited]

Returns:

a reference to the solution of the last weighted sensitivity solve

Definition at line 750 of file system.C.

References System::get_vector().

Referenced by ImplicitSystem::weighted_sensitivity_solve().

{
  return this->get_vector("weighted_sensitivity_solution");
}

bool System::has_variable

(

const std::string &

var

)

const [inherited]

Returns:

true if a variable named var exists in this System

Definition at line 935 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 357 of file implicit_system.h.

References ImplicitSystem::_matrices.

Referenced by ImplicitSystem::add_matrix(), and EigenTimeSolver::init().

{
  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 1444 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, and Threads::spin_mtx.

Referenced by ReferenceCountedObject< SparseMatrix< T > >::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, and Threads::spin_mtx.

Referenced by ReferenceCountedObject< SparseMatrix< T > >::~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 158 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 ImplicitSystem::init_data

(

)

[protected, virtual, inherited]

Initializes the member data fields associated with the system, so that, e.g., assemble() may be used.

Reimplemented from ExplicitSystem.

Reimplemented in ContinuationSystem, DifferentiableSystem, FEMSystem, and FrequencySystem.

Definition at line 85 of file implicit_system.C.

References ImplicitSystem::add_system_matrix(), ExplicitSystem::init_data(), and ImplicitSystem::init_matrices().

Referenced by DifferentiableSystem::init_data().

{
  // initialize parent data
  Parent::init_data();

  // Add the system matrix.
  this->add_system_matrix ();
  
  // Initialize the matrices for the system
  this->init_matrices ();
}

void ImplicitSystem::init_matrices

(

)

[protected, virtual, inherited]

Initializes the matrices associated with this system.

Definition at line 99 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(), 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, \p 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 NewmarkSystem::initial_conditions

(

)

Apply initial conditions.

Definition at line 156 of file newmark_system.C.

References System::get_vector(), and System::user_initialization().

Referenced by assemble().

{
  // libmesh_assert (init_cond_fptr != NULL);

  // Log how long the user's matrix assembly code takes
  START_LOG("initial_conditions ()", "NewmarkSystem");
  
  // Set all values to 0, then
  // call the user-specified function for initial conditions.
  this->get_vector("displacement").zero();
  this->get_vector("velocity").zero();
  this->get_vector("acceleration").zero();
  this->user_initialization();

  // Stop logging the user code
  STOP_LOG("initial_conditions ()", "NewmarkSystem");
}

void System::local_dof_indices	(	const unsigned int	var,
		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 962 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]);
        }
    }
}

unsigned int System::n_active_dofs

(

)

const [inline, inherited]

Returns the number of active degrees of freedom for this System.

Definition at line 1436 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 95 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 88 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 LinearImplicitSystem::n_linear_iterations

(

)

const [inline, inherited]

Returns the number of iterations taken for the most recent linear solve.

Definition at line 140 of file linear_implicit_system.h.

References LinearImplicitSystem::_n_linear_iterations.

{ return _n_linear_iterations; }

unsigned int System::n_local_dofs

(

)

const [inherited]

Returns:

the number of degrees of freedom local to this processor

Definition at line 110 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 364 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.

  { return _n_objects; }

unsigned int System::n_vars

(

)

const [inline, inherited]

Returns:

the number of variables in the system

Definition at line 1390 of file system.h.

References System::_variables.

Referenced by 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(), 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 1452 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 1318 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(), 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 1326 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 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 445 of file system.h.

References System::_solution_projection.

    { return _solution_projection; }

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::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	(	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::prolong_vectors

(

)

[virtual, inherited]

Prolong vectors after the mesh has refined

Definition at line 255 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 ImplicitSystem::qoi_parameter_hessian_vector_product	(	const QoISet &	qoi_indices,
		const ParameterVector &	parameters,
		const ParameterVector &	vector,
		SensitivityData &	product
	)			[virtual, inherited]

For each of the system's quantities of interest q in qoi[qoi_indices], and for a vector of parameters p, the parameter sensitivity Hessian H_ij is defined as H_ij = (d^2 q)/(d p_i d p_j) The Hessian-vector product, for a vector v_k in parameter space, is S_j = H_jk v_k This product is the output of this method, where for each q_i, S_j is stored in sensitivities[i][j].

Reimplemented from System.

Definition at line 825 of file implicit_system.C.

References SensitivityData::allocate_data(), ParameterVector::deep_copy(), QoISet::has_index(), ParameterVector::size(), and ParameterVector::value_copy().

{
  // We currently get partial derivatives via finite differencing
  const Real delta_p = TOLERANCE;

  const unsigned int Np = parameters.size();
  const unsigned int Nq = qoi.size();

  // For each quantity of interest q, the parameter sensitivity
  // Hessian is defined as q''_{kl} = {d^2 q}/{d p_k d p_l}.
  // Given a vector of parameter perturbation weights w_l, this
  // function evaluates the hessian-vector product sum_l(q''_{kl}*w_l)
  //
  // We calculate it from values and partial derivatives of the
  // quantity of interest function Q, solution u, adjoint solution z,
  // parameter sensitivity adjoint solutions z^l, and residual R, as:
  //
  // sum_l(q''_{kl}*w_l) =
  // sum_l(w_l * Q''_{kl}) + Q''_{uk}(u)*(sum_l(w_l u'_l)) -
  // R'_k(u, sum_l(w_l*z^l)) - R'_{uk}(u,z)*(sum_l(w_l u'_l) -
  // sum_l(w_l*R''_{kl}(u,z))
  //
  // See the adjoints model document for more details.

  // We first do an adjoint solve to get z for each quantity of
  // interest

  this->adjoint_solve(qoi_indices);

  // Get ready to fill in senstivities:
  sensitivities.allocate_data(qoi_indices, *this, parameters);

  // We can't solve for all the solution sensitivities u'_l or for all
  // of the parameter sensitivity adjoint solutions z^l without
  // requiring O(Nq*Np) linear solves.  So we'll solve directly for their
  // weighted sum - this is just O(Nq) solves.

  // First solve for sum_l(w_l u'_l).
  this->weighted_sensitivity_solve(parameters, vector);

  // Then solve for sum_l(w_l z^l).
  this->weighted_sensitivity_adjoint_solve(parameters, vector, qoi_indices);

  for (unsigned int k=0; k != Np; ++k)
    {
      // We approximate sum_l(w_l * Q''_{kl}) with a central
      // differencing pertubation:
      // sum_l(w_l * Q''_{kl}) ~= 
      // (Q(p + dp*w_l*e_l + dp*e_k) - Q(p - dp*w_l*e_l + dp*e_k) -
      // Q(p + dp*w_l*e_l - dp*e_k) + Q(p - dp*w_l*e_l - dp*e_k))/(4*dp)

      // The sum(w_l*R''_kl) term requires the same sort of pertubation,
      // and so we subtract it in at the same time:
      // sum_l(w_l * R''_{kl}) ~= 
      // (R(p + dp*w_l*e_l + dp*e_k) - R(p - dp*w_l*e_l + dp*e_k) -
      // R(p + dp*w_l*e_l - dp*e_k) + R(p - dp*w_l*e_l - dp*e_k))/(4*dp)

      ParameterVector oldparameters, parameterperturbation;
      parameters.deep_copy(oldparameters);
      vector.deep_copy(parameterperturbation);
      parameterperturbation *= delta_p;
      parameters += parameterperturbation;

      Number old_parameter = *parameters[k];

      *parameters[k] = old_parameter + delta_p;
      this->assemble_qoi(qoi_indices);
      this->assembly(true, false);
      std::vector<Number> partial2q_term = this->qoi;
      std::vector<Number> partial2R_term(this->qoi.size());
      for (unsigned int i=0; i != Nq; ++i)
        if (qoi_indices.has_index(i))
          partial2R_term[i] = this->rhs->dot(this->get_adjoint_solution(i));

      *parameters[k] = old_parameter - delta_p;
      this->assemble_qoi(qoi_indices);
      this->assembly(true, false);
      for (unsigned int i=0; i != Nq; ++i)
        if (qoi_indices.has_index(i))
          {
            partial2q_term[i] += this->qoi[i];
            partial2R_term[i] += this->rhs->dot(this->get_adjoint_solution(i));
          }

      oldparameters.value_copy(parameters);
      parameterperturbation *= -1.0;
      parameters += parameterperturbation;

      // Re-center old_parameter, which may be affected by vector
      old_parameter = *parameters[k];

      *parameters[k] = old_parameter + delta_p;
      this->assemble_qoi(qoi_indices);
      this->assembly(true, false);
      for (unsigned int i=0; i != Nq; ++i)
        if (qoi_indices.has_index(i))
          {
            partial2q_term[i] += this->qoi[i];
            partial2R_term[i] += this->rhs->dot(this->get_adjoint_solution(i));
          }

      *parameters[k] = old_parameter - delta_p;
      this->assemble_qoi(qoi_indices);
      this->assembly(true, false);
      for (unsigned int i=0; i != Nq; ++i)
        if (qoi_indices.has_index(i))
          {
            partial2q_term[i] += this->qoi[i];
            partial2R_term[i] += this->rhs->dot(this->get_adjoint_solution(i));
          }

      for (unsigned int i=0; i != Nq; ++i)
        if (qoi_indices.has_index(i))
          {
            partial2q_term[i] /= (4. * delta_p);
            partial2R_term[i] /= (4. * delta_p);
          }

      for (unsigned int i=0; i != Nq; ++i)
        if (qoi_indices.has_index(i))
          sensitivities[i][k] = partial2q_term[i];

      // Don't leave the parameters changed
      oldparameters.value_copy(parameters);

      // Re-center old_parameter, which may have been affected by vector
      old_parameter = *parameters[k];

      // We get (partial q / partial u) from the user,
      // but centrally difference it to get q_uk:
      // (partial^2 q / partial u partial k)
      // q_uk = (q_u(p+dp*e_k) - q_u(p-dp*e_k))/(2*dp)

      // To avoid creating Nq temporary vectors, we add these terms to
      // the sensitivities output one by one.
      //
      // FIXME: this is probably a bad order of operations for
      // controlling floating point error.

      *parameters[k] = old_parameter + delta_p;
      this->assemble_qoi_derivative(qoi_indices);

      for (unsigned int i=0; i != Nq; ++i)
        if (qoi_indices.has_index(i))
          sensitivities[i][k] += this->get_adjoint_rhs(i).dot(this->get_weighted_sensitivity_solution()) / (2.*delta_p);
 
      *parameters[k] = old_parameter - delta_p;
      this->assemble_qoi_derivative(qoi_indices);

      for (unsigned int i=0; i != Nq; ++i)
        if (qoi_indices.has_index(i))
          sensitivities[i][k] -= this->get_adjoint_rhs(i).dot(this->get_weighted_sensitivity_solution()) / (2.*delta_p);

      // We get R from the user,
      // but centrally difference it to get R_k:
      // (partial R / partial k)
      // R_k*sum(w_l*z^l) = (R(p+dp*e_k)*sum(w_l*z^l) - R(p-dp*e_k)*sum(w_l*z^l))/(2*dp)
      *parameters[k] = old_parameter + delta_p;
      this->assembly(true, false);

      std::vector<Number> partialR_term(this->qoi.size());
      for (unsigned int i=0; i != Nq; ++i)
        if (qoi_indices.has_index(i))
          partialR_term[i] = this->rhs->dot(this->get_weighted_sensitivity_adjoint_solution(i));

      *parameters[k] = old_parameter - delta_p;
      this->assembly(true, false);

      for (unsigned int i=0; i != Nq; ++i)
        if (qoi_indices.has_index(i))
          {
            partialR_term[i] -= this->rhs->dot(this->get_weighted_sensitivity_adjoint_solution(i));
            partialR_term[i] /= (2.*delta_p);

            // We subtract since this term is negative
            sensitivities[i][k] -= partialR_term[i];
          }

      // We get (partial R / partial u) from the user,
      // but centrally difference it to get R_uk*z*sum(w_l*u'_l):
      // (partial^2 R / partial u partial k)
      // R_uk*z = (R_u(p+dp*e_k)*z*sum(w_l*u'_l) - R_u(p-dp*e_k)*z*sum(w_l*u'_l))/(2*dp)

      // To avoid creating Nq temporary vectors, we again add these
      // terms to the sensitivities output one by one.
      //
      // FIXME: this is probably a bad order of operations for
      // controlling floating point error.

      *parameters[k] = old_parameter + delta_p;
      this->assembly(false, true);

      // Use a single temporary vector for matrix-vector-vector products
      AutoPtr<NumericVector<Number> > tempvec = this->get_weighted_sensitivity_solution().zero_clone();
      this->matrix->vector_mult(*tempvec, this->get_weighted_sensitivity_solution());
      for (unsigned int i=0; i != Nq; ++i)
        if (qoi_indices.has_index(i))
          {
            // We subtract for the +dp component since this term is negative
            sensitivities[i][k] -= 
              this->get_adjoint_solution(i).dot(*tempvec) / (2.*delta_p);
          }

      *parameters[k] = old_parameter - delta_p;
      this->assembly(false, true);

      this->matrix->vector_mult(*tempvec, this->get_weighted_sensitivity_solution());
      for (unsigned int i=0; i != Nq; ++i)
        if (qoi_indices.has_index(i))
          {
            sensitivities[i][k] += 
              this->get_adjoint_solution(i).dot(*tempvec) / (2.*delta_p);
          }
    }

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

void System::qoi_parameter_sensitivity	(	const QoISet &	qoi_indices,
		const ParameterVector &	parameters,
		SensitivityData &	sensitivities
	)			[virtual, inherited]

Solves for the derivative of each of the system's quantities of interest q in qoi[qoi_indices] with respect to each parameter in parameters, placing the result for qoi i and parameter j into sensitivities[i][j].

Note that parameters is a const vector, not a vector-of-const; parameter values in this vector need to be mutable for finite differencing to work.

Automatically chooses the forward method for problems with more quantities of interest than parameters, or the adjoint method otherwise.

This method is only usable in derived classes which overload an implementation.

Definition at line 383 of file system.C.

References ParameterVector::size(), and QoISet::size().

{
  // Forward sensitivities are more efficient for Nq > Np
  if (qoi_indices.size(*this) > parameters.size())
    forward_qoi_parameter_sensitivity(qoi_indices, parameters, sensitivities);
  // Adjoint sensitivities are more efficient for Np > Nq,
  // and an adjoint may be more reusable than a forward 
  // solution sensitivity in the Np == Nq case.
  else
    adjoint_qoi_parameter_sensitivity(qoi_indices, parameters, sensitivities);
}

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 314 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 bool	read_header = true,
		const bool	read_additional_data = true,
		const bool	read_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(), type, 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 << "*****************************************************************\n"
                        << "* WARNING: reading a potentially incompatible restart file!!!   *\n"
                        << "*  contact libmesh-users@lists.sourceforge.net for more details *\n"
                        << "*****************************************************************"
                        << 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 bool	read_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 bool	read_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 bool	read_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 NewmarkSystem::reinit

(

)

[virtual]

Reinitializes the member data fields associated with the system, so that, e.g., assemble() may be used.

Reimplemented from LinearImplicitSystem.

Definition at line 127 of file newmark_system.C.

{
  libmesh_error();
  
  // initialize parent data
  LinearImplicitSystem::reinit();
}

void LinearImplicitSystem::release_linear_solver

(

LinearSolver< Number > *

s

)

const [virtual, inherited]

Releases a pointer to a linear solver acquired by this->get_linear_solver()

Reimplemented from ImplicitSystem.

Definition at line 327 of file linear_implicit_system.C.

{
}

SparseMatrix< Number > * ImplicitSystem::request_matrix

(

const std::string &

mat_name

)

[inherited]

Returns:

a writable pointer to this system's additional matrix named mat_name, or returns NULL if no matrix by that name exists.

Definition at line 213 of file implicit_system.C.

References ImplicitSystem::_matrices.

{
  // Make sure the matrix exists
  matrices_iterator pos = _matrices.find (mat_name);
  
  if (pos == _matrices.end())
    return NULL;
  
  return pos->second;
}

const SparseMatrix< Number > * ImplicitSystem::request_matrix

(

const std::string &

mat_name

)

const [inherited]

Returns:

a const pointer to this system's additional matrix named mat_name, or returns NULL if no matrix by that name exists.

Definition at line 200 of file implicit_system.C.

References ImplicitSystem::_matrices.

Referenced by ImplicitSystem::sensitivity_solve(), NewtonSolver::solve(), and LinearImplicitSystem::solve().

{
  // Make sure the matrix exists
  const_matrices_iterator pos = _matrices.find (mat_name);
  
  if (pos == _matrices.end())
    return NULL;
  
  return pos->second;
}

NumericVector< Number > * System::request_vector

(

const unsigned int

vec_num

)

[inherited]

Returns:

a writeable pointer to this system's additional vector number vec_num (where the vectors are counted starting with 0), or returns NULL if the system has no such vector.

Definition at line 611 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;
    }
  if (v==v_end)
    return NULL;
  return v->second;
}

const NumericVector< Number > * System::request_vector

(

const unsigned int

vec_num

)

const [inherited]

Returns:

a const pointer to this system's additional vector number vec_num (where the vectors are counted starting with 0), or returns NULL if the system has no such vector.

Definition at line 594 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;
    }
  if (v==v_end)
    return NULL;
  return v->second;
}

NumericVector< Number > * System::request_vector

(

const std::string &

vec_name

)

[inherited]

Returns:

a pointer to the vector if this System has a vector associated with the given name, NULL otherwise.

Definition at line 582 of file system.C.

References System::_vectors.

{
  vectors_iterator pos = _vectors.find(vec_name);
  
  if (pos == _vectors.end())
    return NULL;
  
  return pos->second;
}

const NumericVector< Number > * System::request_vector

(

const std::string &

vec_name

)

const [inherited]

Returns:

a const pointer to the vector if this System has a vector associated with the given name, NULL otherwise.

Definition at line 570 of file system.C.

References System::_vectors.

{
  const_vectors_iterator pos = _vectors.find(vec_name);
  
  if (pos == _vectors.end())
    return NULL;
  
  return pos->second;
}

void System::restrict_vectors

(

)

[virtual, inherited]

Restrict vectors after the mesh has coarsened

Definition at line 219 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
}

std::pair< unsigned int, Real > ImplicitSystem::sensitivity_solve

(

const ParameterVector &

parameters

)

[virtual, inherited]

Assembles & solves the linear system(s) (dR/du)*u_p = -dR/dp, for those parameters contained within parameters.

Returns a pair with the total number of linear iterations performed and the (sum of the) final residual norms

Reimplemented from System.

Definition at line 282 of file implicit_system.C.

References System::assemble_before_solve, ImplicitSystem::assemble_residual_derivatives(), ImplicitSystem::assembly(), DofMap::enforce_constraints_exactly(), System::get_dof_map(), ImplicitSystem::get_linear_solve_parameters(), ImplicitSystem::get_linear_solver(), System::get_sensitivity_rhs(), System::get_sensitivity_solution(), ImplicitSystem::matrix, ImplicitSystem::release_linear_solver(), ImplicitSystem::request_matrix(), ParameterVector::size(), and LinearSolver< T >::solve().

{
  // Log how long the linear solve takes.
  START_LOG("sensitivity_solve()", "ImplicitSystem");

  // The forward system should now already be solved.
  // Now assemble the corresponding sensitivity system.

  if (this->assemble_before_solve)
    {
      // Build the Jacobian
      this->assembly(false, true);
      this->matrix->close();

      // Reset and build the RHS from the residual derivatives
      this->assemble_residual_derivatives(parameters);
    }

  // The sensitivity problem is linear
  LinearSolver<Number> *linear_solver = this->get_linear_solver();

  // Our iteration counts and residuals will be sums of the individual
  // results
  std::pair<unsigned int, Real> solver_params =
    this->get_linear_solve_parameters();
  std::pair<unsigned int, Real> totalrval = std::make_pair(0,0.0);

  // Solve the linear system.
  SparseMatrix<Number> *pc = this->