libMesh::LinearImplicitSystem Class Reference

#include <linear_implicit_system.h>

Inheritance diagram for libMesh::LinearImplicitSystem:

Public Types
typedef LinearImplicitSystem	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
	LinearImplicitSystem (EquationSystems &es, const std::string &name, const unsigned int number)
virtual	~LinearImplicitSystem ()
sys_type &	system ()
virtual void	clear ()
virtual void	init_data ()
virtual void	reinit ()
virtual void	assemble ()
virtual void	restrict_solve_to (const SystemSubset *subset, const SubsetSolveMode subset_solve_mode=SUBSET_ZERO)
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)
virtual std::string	system_type () const
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 void	disable_cache ()
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 (const QoISet &qoi_indices, const ParameterVector &parameters, SensitivityData &hessian)
virtual void	qoi_parameter_hessian_vector_product (const QoISet &qoi_indices, const ParameterVector &parameters, const ParameterVector &vector, SensitivityData &product)
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)
virtual 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 ()
bool	is_adjoint_already_solved () const
void	set_adjoint_already_solved (bool setting)
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 (FunctionBase< Number > f, FunctionBase< Gradient > g=NULL) const
void	project_solution (FEMFunctionBase< Number > f, FEMFunctionBase< Gradient > g=NULL) const
void	project_solution (Number fptr(const Point &p, const Parameters &parameters, const std::string &sys_name, const std::string &unknown_name), Gradient gptr(const Point &p, const Parameters &parameters, const std::string &sys_name, const std::string &unknown_name), const Parameters &parameters) const
void	project_vector (NumericVector< Number > &new_vector, FunctionBase< Number > f, FunctionBase< Gradient > g=NULL) const
void	project_vector (NumericVector< Number > &new_vector, FEMFunctionBase< Number > f, FEMFunctionBase< Gradient > g=NULL) const
void	project_vector (Number fptr(const Point &p, const Parameters &parameters, const std::string &sys_name, const std::string &unknown_name), Gradient gptr(const Point &p, const Parameters &parameters, const std::string &sys_name, const std::string &unknown_name), const Parameters &parameters, NumericVector< Number > &new_vector) const
void	boundary_project_solution (const std::set< boundary_id_type > &b, const std::vector< unsigned int > &variables, FunctionBase< Number > f, FunctionBase< Gradient > g=NULL)
void	boundary_project_solution (const std::set< boundary_id_type > &b, const std::vector< unsigned int > &variables, Number fptr(const Point &p, const Parameters &parameters, const std::string &sys_name, const std::string &unknown_name), Gradient gptr(const Point &p, const Parameters &parameters, const std::string &sys_name, const std::string &unknown_name), const Parameters &parameters)
void	boundary_project_vector (const std::set< boundary_id_type > &b, const std::vector< unsigned int > &variables, NumericVector< Number > &new_vector, FunctionBase< Number > f, FunctionBase< Gradient > g=NULL) const
void	boundary_project_vector (const std::set< boundary_id_type > &b, const std::vector< unsigned int > &variables, Number fptr(const Point &p, const Parameters &parameters, const std::string &sys_name, const std::string &unknown_name), Gradient gptr(const Point &p, const Parameters &parameters, const std::string &sys_name, const std::string &unknown_name), const Parameters &parameters, NumericVector< Number > &new_vector) const
unsigned int	number () const
void	update_global_solution (std::vector< Number > &global_soln) const
void	update_global_solution (std::vector< Number > &global_soln, const unsigned int dest_proc) const
const MeshBase &	get_mesh () const
MeshBase &	get_mesh ()
const DofMap &	get_dof_map () const
DofMap &	get_dof_map ()
const EquationSystems &	get_equation_systems () const
EquationSystems &	get_equation_systems ()
bool	active () const
void	activate ()
void	deactivate ()
void	set_basic_system_only ()
vectors_iterator	vectors_begin ()
const_vectors_iterator	vectors_begin () const
vectors_iterator	vectors_end ()
const_vectors_iterator	vectors_end () const
NumericVector< Number > &	add_vector (const std::string &vec_name, const bool projections=true, const ParallelType type=PARALLEL)
void	remove_vector (const std::string &vec_name)
bool &	project_solution_on_reinit (void)
bool	have_vector (const std::string &vec_name) const
const NumericVector< Number > *	request_vector (const std::string &vec_name) const
NumericVector< Number > *	request_vector (const std::string &vec_name)
const NumericVector< Number > *	request_vector (const unsigned int vec_num) const
NumericVector< Number > *	request_vector (const unsigned int vec_num)
const NumericVector< Number > &	get_vector (const std::string &vec_name) const
NumericVector< Number > &	get_vector (const std::string &vec_name)
const NumericVector< Number > &	get_vector (const unsigned int vec_num) const
NumericVector< Number > &	get_vector (const unsigned int vec_num)
const std::string &	vector_name (const unsigned int vec_num) const
const std::string &	vector_name (const NumericVector< Number > &vec_reference) const
void	set_vector_preservation (const std::string &vec_name, bool preserve)
bool	vector_preservation (const std::string &vec_name) const
NumericVector< Number > &	add_adjoint_solution (unsigned int i=0)
NumericVector< Number > &	get_adjoint_solution (unsigned int i=0)
const NumericVector< Number > &	get_adjoint_solution (unsigned int i=0) const
NumericVector< Number > &	add_sensitivity_solution (unsigned int i=0)
NumericVector< Number > &	get_sensitivity_solution (unsigned int i=0)
const NumericVector< Number > &	get_sensitivity_solution (unsigned int i=0) const
NumericVector< Number > &	add_weighted_sensitivity_adjoint_solution (unsigned int i=0)
NumericVector< Number > &	get_weighted_sensitivity_adjoint_solution (unsigned int i=0)
const NumericVector< Number > &	get_weighted_sensitivity_adjoint_solution (unsigned int i=0) const
NumericVector< Number > &	add_weighted_sensitivity_solution ()
NumericVector< Number > &	get_weighted_sensitivity_solution ()
const NumericVector< Number > &	get_weighted_sensitivity_solution () const
NumericVector< Number > &	add_adjoint_rhs (unsigned int i=0)
NumericVector< Number > &	get_adjoint_rhs (unsigned int i=0)
const NumericVector< Number > &	get_adjoint_rhs (unsigned int i=0) const
NumericVector< Number > &	add_sensitivity_rhs (unsigned int i=0)
NumericVector< Number > &	get_sensitivity_rhs (unsigned int i=0)
const NumericVector< Number > &	get_sensitivity_rhs (unsigned int i=0) const
unsigned int	n_vectors () const
unsigned int	n_vars () const
unsigned int	n_variable_groups () const
unsigned int	n_components () const
dof_id_type	n_dofs () const
dof_id_type	n_active_dofs () const
dof_id_type	n_constrained_dofs () const
dof_id_type	n_local_constrained_dofs () const
dof_id_type	n_local_dofs () const
unsigned int	add_variable (const std::string &var, const FEType &type, const std::set< subdomain_id_type > *const active_subdomains=NULL)
unsigned int	add_variable (const std::string &var, const Order order=FIRST, const FEFamily=LAGRANGE, const std::set< subdomain_id_type > *const active_subdomains=NULL)
unsigned int	add_variables (const std::vector< std::string > &vars, const FEType &type, const std::set< subdomain_id_type > *const active_subdomains=NULL)
unsigned int	add_variables (const std::vector< std::string > &vars, const Order order=FIRST, const FEFamily=LAGRANGE, const std::set< subdomain_id_type > *const active_subdomains=NULL)
const Variable &	variable (unsigned int var) const
const VariableGroup &	variable_group (unsigned int vg) const
bool	has_variable (const std::string &var) const
const std::string &	variable_name (const unsigned int i) const
unsigned short int	variable_number (const std::string &var) const
void	get_all_variable_numbers (std::vector< unsigned int > &all_variable_numbers) const
unsigned int	variable_scalar_number (const std::string &var, unsigned int component) const
unsigned int	variable_scalar_number (unsigned int var_num, unsigned int component) const
const FEType &	variable_type (const unsigned int i) const
const FEType &	variable_type (const std::string &var) const
bool	identify_variable_groups () const
void	identify_variable_groups (const bool)
Real	calculate_norm (const NumericVector< Number > &v, unsigned int var=0, FEMNormType norm_type=L2) const
Real	calculate_norm (const NumericVector< Number > &v, const SystemNorm &norm) const
void	read_header (Xdr &io, const std::string &version, const bool read_header=true, const bool read_additional_data=true, const bool read_legacy_format=false)
void	read_legacy_data (Xdr &io, const bool read_additional_data=true)
void	read_serialized_data (Xdr &io, const bool read_additional_data=true)
dof_id_type	read_serialized_vectors (Xdr &io, const std::vector< NumericVector< Number > * > &vectors) const
void	read_parallel_data (Xdr &io, const bool read_additional_data)
void	write_header (Xdr &io, const std::string &version, const bool write_additional_data) const
void	write_serialized_data (Xdr &io, const bool write_additional_data=true) const
dof_id_type	write_serialized_vectors (Xdr &io, const std::vector< const NumericVector< Number > * > &vectors) const
void	write_parallel_data (Xdr &io, const bool write_additional_data) const
std::string	get_info () const
void	attach_init_function (void fptr(EquationSystems &es, const std::string &name))
void	attach_init_object (Initialization &init)
void	attach_assemble_function (void fptr(EquationSystems &es, const std::string &name))
void	attach_assemble_object (Assembly &assemble)
void	attach_constraint_function (void fptr(EquationSystems &es, const std::string &name))
void	attach_constraint_object (Constraint &constrain)
void	attach_QOI_function (void fptr(EquationSystems &es, const std::string &name, const QoISet &qoi_indices))
void	attach_QOI_object (QOI &qoi)
void	attach_QOI_derivative (void fptr(EquationSystems &es, const std::string &name, const QoISet &qoi_indices))
void	attach_QOI_derivative_object (QOIDerivative &qoi_derivative)
virtual void	user_initialization ()
virtual void	user_assembly ()
virtual void	user_constrain ()
virtual void	user_QOI (const QoISet &qoi_indices)
virtual void	user_QOI_derivative (const QoISet &qoi_indices)
virtual void	re_update ()
virtual void	restrict_vectors ()
virtual void	prolong_vectors ()
Number	current_solution (const dof_id_type global_dof_number) const
Number	point_value (unsigned int var, const Point &p, const bool insist_on_success=true) const
Number	point_value (unsigned int var, const Point &p, const Elem &e) const
Gradient	point_gradient (unsigned int var, const Point &p, const bool insist_on_success=true) const
Gradient	point_gradient (unsigned int var, const Point &p, const Elem &e) const
Tensor	point_hessian (unsigned int var, const Point &p, const bool insist_on_success=true) const
Tensor	point_hessian (unsigned int var, const Point &p, const Elem &e) const
void	local_dof_indices (const unsigned int var, std::set< dof_id_type > &var_indices) const
void	zero_variable (NumericVector< Number > &v, unsigned int var_num) const
Static Public Member Functions
static std::string	get_info ()
static void	print_info (std::ostream &out=libMesh::out)
static unsigned int	n_objects ()
static void	enable_print_counter_info ()
static void	disable_print_counter_info ()
Public Attributes
AutoPtr< LinearSolver< Number > >	linear_solver
SparseMatrix< Number > *	matrix
NumericVector< Number > *	rhs
bool	assemble_before_solve
bool	use_fixed_solution
int	extra_quadrature_order
AutoPtr< NumericVector< Number > >	solution
AutoPtr< NumericVector< Number > >	current_local_solution
Real	time
std::vector< Number >	qoi
Protected Types
typedef std::map< std::string, std::pair< unsigned int, unsigned int > >	Counts
Protected Member Functions
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
const SystemSubset *	_subset
SubsetSolveMode	_subset_solve_mode
Static Protected Attributes
static Counts	_counts
static Threads::atomic < unsigned int >	_n_objects
static Threads::spin_mutex	_mutex
static bool	_enable_print_counter = true

Detailed Description

This class provides a specific system class. It aims at implicit systems, offering nothing more than just the essentials needed to solve a system. Note that still additional vectors/matrices may be added, as offered in the parent class ExplicitSystem.

Definition at line 49 of file linear_implicit_system.h.

Member Typedef Documentation

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

Definition at line 277 of file implicit_system.h.

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

Definition at line 716 of file system.h.

typedef std::map<std::string, std::pair<unsigned int, unsigned int> > libMesh::ReferenceCounter::Counts [protected, inherited]

Data structure to log the information. The log is identified by the class name.

Definition at line 113 of file reference_counter.h.

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

Matrix iterator typedefs.

Definition at line 276 of file implicit_system.h.

typedef ImplicitSystem libMesh::LinearImplicitSystem::Parent

The type of the parent.

Reimplemented from libMesh::ImplicitSystem.

Reimplemented in libMesh::RBConstruction, libMesh::RBEIMConstruction, and libMesh::TransientRBConstruction.

Definition at line 74 of file linear_implicit_system.h.

typedef LinearImplicitSystem libMesh::LinearImplicitSystem::sys_type

The type of system.

Reimplemented from libMesh::ImplicitSystem.

Reimplemented in libMesh::NewmarkSystem, libMesh::RBConstruction, libMesh::RBEIMConstruction, libMesh::TransientRBConstruction, libMesh::RBConstructionBase< LinearImplicitSystem >, and libMesh::TransientSystem< RBConstruction >.

Definition at line 69 of file linear_implicit_system.h.

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

Vector iterator typedefs.

Definition at line 715 of file system.h.

Constructor & Destructor Documentation

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

Constructor. Optionally initializes required data structures.

Definition at line 37 of file linear_implicit_system.C.

00039                                                                           :
00040 
00041   Parent                 (es, name_in, number_in),
00042   linear_solver          (LinearSolver<Number>::build()),
00043   _n_linear_iterations   (0),
00044   _final_linear_residual (1.e20),
00045   _shell_matrix(NULL),
00046   _subset(NULL),
00047   _subset_solve_mode(SUBSET_ZERO)
00048 {
00049 }

libMesh::LinearImplicitSystem::~LinearImplicitSystem ( ) [virtual]

Destructor.

Definition at line 53 of file linear_implicit_system.C.

References clear().

00054 {
00055   // Clear data
00056   this->clear();
00057 }

Member Function Documentation

void libMesh::System::activate ( ) [inline, inherited]

Activates the system. Only active systems are solved.

Definition at line 1874 of file system.h.

References libMesh::System::_active.

01875 {
01876   _active = true;
01877 }

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

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

Definition at line 1866 of file system.h.

References libMesh::System::_active.

01867 {
01868   return _active;
01869 }

NumericVector< Number > & libMesh::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 1016 of file system.C.

References libMesh::System::add_vector().

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

01017 {
01018   std::ostringstream adjoint_rhs_name;
01019   adjoint_rhs_name << "adjoint_rhs" << i;
01020 
01021   return this->add_vector(adjoint_rhs_name.str(), false);
01022 }

NumericVector< Number > & libMesh::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 956 of file system.C.

References libMesh::System::add_vector().

Referenced by libMesh::ImplicitSystem::adjoint_solve().

00957 {
00958   std::ostringstream adjoint_name;
00959   adjoint_name << "adjoint_solution" << i;
00960 
00961   return this->add_vector(adjoint_name.str());
00962 }

SparseMatrix< Number > & libMesh::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 207 of file implicit_system.C.

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

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

00208 {
00209   // only add matrices before initializing...
00210   if (!_can_add_matrices)
00211     {
00212       libMesh::err << "ERROR: Too late.  Cannot add matrices to the system after initialization"
00213                     << std::endl
00214                     << " any more.  You should have done this earlier."
00215                     << std::endl;
00216       libmesh_error();
00217     }
00218 
00219   // Return the matrix if it is already there.
00220   if (this->have_matrix(mat_name))
00221     return *(_matrices[mat_name]);
00222 
00223   // Otherwise build the matrix and return it.
00224   SparseMatrix<Number>* buf = SparseMatrix<Number>::build().release();
00225   _matrices.insert (std::make_pair (mat_name, buf));
00226 
00227   return *buf;
00228 }

NumericVector< Number > & libMesh::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 1046 of file system.C.

References libMesh::System::add_vector().

Referenced by libMesh::ImplicitSystem::assemble_residual_derivatives().

01047 {
01048   std::ostringstream sensitivity_rhs_name;
01049   sensitivity_rhs_name << "sensitivity_rhs" << i;
01050 
01051   return this->add_vector(sensitivity_rhs_name.str(), false);
01052 }

NumericVector< Number > & libMesh::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 905 of file system.C.

References libMesh::System::add_vector().

Referenced by libMesh::ImplicitSystem::sensitivity_solve().

00906 {
00907   std::ostringstream sensitivity_name;
00908   sensitivity_name << "sensitivity_solution" << i;
00909 
00910   return this->add_vector(sensitivity_name.str());
00911 }

unsigned int libMesh::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 1155 of file system.C.

References libMesh::System::add_variable().

01159 {
01160   return this->add_variable(var,
01161                             FEType(order, family),
01162                             active_subdomains);
01163 }

unsigned int libMesh::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 1076 of file system.C.

References libMesh::System::_variable_groups, libMesh::System::_variable_numbers, libMesh::System::_variables, libMesh::System::add_variables(), libMesh::err, libMesh::System::identify_variable_groups(), libMesh::System::n_variable_groups(), libMesh::System::n_vars(), libMesh::System::number(), libMesh::System::variable_name(), and libMesh::System::variable_type().

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

01079 {
01080   // Make sure the variable isn't there already
01081   // or if it is, that it's the type we want
01082   for (unsigned int v=0; v<this->n_vars(); v++)
01083     if (this->variable_name(v) == var)
01084       {
01085         if (this->variable_type(v) == type)
01086           return _variables[v].number();
01087 
01088         libMesh::err << "ERROR: incompatible variable "
01089                      << var
01090                      << " has already been added for this system!"
01091                      << std::endl;
01092         libmesh_error();
01093       }
01094 
01095   // Optimize for VariableGroups here - if the user is adding multiple
01096   // variables of the same FEType and subdomain restriction, catch
01097   // that here and add them as members of the same VariableGroup.
01098   //
01099   // start by setting this flag to whatever the user has requested
01100   // and then consider the conditions which should negate it.
01101   bool should_be_in_vg = this->identify_variable_groups();
01102 
01103   // No variable groups, nothing to add to
01104   if (!this->n_variable_groups())
01105     should_be_in_vg = false;
01106 
01107   else
01108     {
01109       VariableGroup &vg(_variable_groups.back());
01110 
01111       // get a pointer to their subdomain restriction, if any.
01112       const std::set<subdomain_id_type> * const
01113         their_active_subdomains (vg.implicitly_active() ?
01114                                  NULL : &vg.active_subdomains());
01115 
01116       // Different types?
01117       if (vg.type() != type)
01118         should_be_in_vg = false;
01119 
01120       // they are restricted, we aren't?
01121       if (their_active_subdomains && !active_subdomains)
01122         should_be_in_vg = false;
01123 
01124       // they aren't restriced, we are?
01125       if (!their_active_subdomains && active_subdomains)
01126         should_be_in_vg = false;
01127 
01128       if (their_active_subdomains && active_subdomains)
01129         // restricted to different sets?
01130         if (*their_active_subdomains != *active_subdomains)
01131           should_be_in_vg = false;
01132 
01133       // OK, after all that, append the variable to the vg if none of the conditions
01134       // were violated
01135       if (should_be_in_vg)
01136         {
01137           const unsigned int curr_n_vars = this->n_vars();
01138 
01139           vg.append (var);
01140 
01141           _variables.push_back(vg(vg.n_variables()-1));
01142           _variable_numbers[var] = curr_n_vars;
01143           return curr_n_vars;
01144         }
01145     }
01146 
01147   // otherwise, fall back to adding a single variable group
01148   return this->add_variables (std::vector<std::string>(1, var),
01149                               type,
01150                               active_subdomains);
01151 }

unsigned int libMesh::System::add_variables	(	const std::vector< std::string > &	vars,
		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 1220 of file system.C.

References libMesh::System::add_variables().

01224 {
01225   return this->add_variables(vars,
01226                              FEType(order, family),
01227                              active_subdomains);
01228 }

unsigned int libMesh::System::add_variables	(	const std::vector< std::string > &	vars,
		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 1167 of file system.C.

References libMesh::System::_variable_groups, libMesh::System::_variable_numbers, libMesh::System::_variables, libMesh::err, libMesh::System::n_components(), libMesh::System::n_vars(), libMesh::System::number(), libMesh::System::variable_name(), and libMesh::System::variable_type().

Referenced by libMesh::System::add_variable(), and libMesh::System::add_variables().

01170 {
01171   // Make sure the variable isn't there already
01172   // or if it is, that it's the type we want
01173   for (unsigned int ov=0; ov<vars.size(); ov++)
01174     for (unsigned int v=0; v<this->n_vars(); v++)
01175       if (this->variable_name(v) == vars[ov])
01176         {
01177           if (this->variable_type(v) == type)
01178             return _variables[v].number();
01179 
01180           libMesh::err << "ERROR: incompatible variable "
01181                        << vars[ov]
01182                        << " has already been added for this system!"
01183                        << std::endl;
01184           libmesh_error();
01185         }
01186 
01187   const unsigned int curr_n_vars = this->n_vars();
01188 
01189   const unsigned int next_first_component = this->n_components();
01190 
01191   // Add the variable group to the list
01192   _variable_groups.push_back((active_subdomains == NULL) ?
01193                              VariableGroup(this, vars, curr_n_vars,
01194                                            next_first_component, type) :
01195                              VariableGroup(this, vars, curr_n_vars,
01196                                            next_first_component, type, *active_subdomains));
01197 
01198   const VariableGroup &vg (_variable_groups.back());
01199 
01200   // Add each component of the group individually
01201   for (unsigned int v=0; v<vars.size(); v++)
01202     {
01203       _variables.push_back (vg(v));
01204       _variable_numbers[vars[v]] = curr_n_vars+v;
01205     }
01206 
01207   libmesh_assert_equal_to ((curr_n_vars+vars.size()), this->n_vars());
01208 
01209   // BSK - Defer this now to System::init_data() so we can detect
01210   // VariableGroups 12/28/2012
01211   // // Add the variable group to the _dof_map
01212   // _dof_map->add_variable_group (vg);
01213 
01214   // Return the number of the new variable
01215   return curr_n_vars+vars.size()-1;
01216 }

NumericVector< Number > & libMesh::System::add_vector	(	const std::string &	vec_name,
		const bool	projections = `true`,
		const ParallelType	type = `PARALLEL`
	)			`[inherited]`

Adds the additional vector vec_name to this system. 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 675 of file system.C.

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

Referenced by libMesh::System::add_adjoint_rhs(), libMesh::System::add_adjoint_solution(), libMesh::System::add_sensitivity_rhs(), libMesh::System::add_sensitivity_solution(), libMesh::ExplicitSystem::add_system_rhs(), libMesh::System::add_weighted_sensitivity_adjoint_solution(), libMesh::System::add_weighted_sensitivity_solution(), libMesh::UnsteadySolver::init(), libMesh::ContinuationSystem::init_data(), libMesh::NewmarkSystem::NewmarkSystem(), libMesh::System::read_header(), libMesh::FrequencySystem::set_frequencies(), libMesh::FrequencySystem::set_frequencies_by_range(), and libMesh::FrequencySystem::set_frequencies_by_steps().

00678 {
00679   // Return the vector if it is already there.
00680   if (this->have_vector(vec_name))
00681     return *(_vectors[vec_name]);
00682 
00683   // Otherwise build the vector
00684   NumericVector<Number>* buf = NumericVector<Number>::build().release();
00685   _vectors.insert (std::make_pair (vec_name, buf));
00686   _vector_projections.insert (std::make_pair (vec_name, projections));
00687 
00688   _vector_types.insert (std::make_pair (vec_name, type));
00689 
00690   // Initialize it if necessary
00691   if (!_can_add_vectors)
00692   {
00693     if(type == GHOSTED)
00694     {
00695 #ifdef LIBMESH_ENABLE_GHOSTED
00696       buf->init (this->n_dofs(), this->n_local_dofs(),
00697                  _dof_map->get_send_list(), false,
00698                  GHOSTED);
00699 #else
00700       std::cerr << "Cannot initialize ghosted vectors when they are not enabled." << std::endl;
00701       libmesh_error();
00702 #endif
00703     }
00704     else
00705       buf->init (this->n_dofs(), this->n_local_dofs(), false, type);
00706   }
00707 
00708   return *buf;
00709 }

NumericVector< Number > & libMesh::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 986 of file system.C.

References libMesh::System::add_vector().

Referenced by libMesh::ImplicitSystem::weighted_sensitivity_adjoint_solve().

00987 {
00988   std::ostringstream adjoint_name;
00989   adjoint_name << "weighted_sensitivity_adjoint_solution" << i;
00990 
00991   return this->add_vector(adjoint_name.str());
00992 }

NumericVector< Number > & libMesh::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 935 of file system.C.

References libMesh::System::add_vector().

Referenced by libMesh::ImplicitSystem::weighted_sensitivity_solve().

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

void libMesh::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 libMesh::System.

Definition at line 696 of file implicit_system.C.

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

00699 {
00700   const unsigned int Np = libmesh_cast_int<unsigned int>
00701     (parameters.size());
00702   const unsigned int Nq = libmesh_cast_int<unsigned int>
00703     (qoi.size());
00704 
00705   // We currently get partial derivatives via central differencing
00706   const Real delta_p = TOLERANCE;
00707 
00708   // An introduction to the problem:
00709   //
00710   // Residual R(u(p),p) = 0
00711   // partial R / partial u = J = system matrix
00712   //
00713   // This implies that:
00714   // d/dp(R) = 0
00715   // (partial R / partial p) +
00716   // (partial R / partial u) * (partial u / partial p) = 0
00717 
00718   // We first do an adjoint solve:
00719   // J^T * z = (partial q / partial u)
00720   // if we havent already or dont have an initial condition for the adjoint
00721   if (!this->is_adjoint_already_solved())
00722     {
00723       this->adjoint_solve(qoi_indices);
00724     }
00725 
00726   // Get ready to fill in senstivities:
00727   sensitivities.allocate_data(qoi_indices, *this, parameters);
00728 
00729   // We use the identities:
00730   // dq/dp = (partial q / partial p) + (partial q / partial u) *
00731   //         (partial u / partial p)
00732   // dq/dp = (partial q / partial p) + (J^T * z) *
00733   //         (partial u / partial p)
00734   // dq/dp = (partial q / partial p) + z * J *
00735   //         (partial u / partial p)
00736 
00737   // Leading to our final formula:
00738   // dq/dp = (partial q / partial p) - z * (partial R / partial p)
00739 
00740   for (unsigned int j=0; j != Np; ++j)
00741     {
00742       // (partial q / partial p) ~= (q(p+dp)-q(p-dp))/(2*dp)
00743       // (partial R / partial p) ~= (rhs(p+dp) - rhs(p-dp))/(2*dp)
00744 
00745       Number old_parameter = *parameters[j];
00746       // Number old_qoi = this->qoi;
00747 
00748       *parameters[j] = old_parameter - delta_p;
00749       this->assemble_qoi(qoi_indices);
00750       std::vector<Number> qoi_minus = this->qoi;
00751 
00752       this->assembly(true, false);
00753       this->rhs->close();
00754 
00755       // FIXME - this can and should be optimized to avoid the clone()
00756       AutoPtr<NumericVector<Number> > partialR_partialp = this->rhs->clone();
00757       *partialR_partialp *= -1;
00758 
00759       *parameters[j] = old_parameter + delta_p;
00760       this->assemble_qoi(qoi_indices);
00761       std::vector<Number>& qoi_plus = this->qoi;
00762 
00763       std::vector<Number> partialq_partialp(Nq, 0);
00764       for (unsigned int i=0; i != Nq; ++i)
00765         if (qoi_indices.has_index(i))
00766           partialq_partialp[i] = (qoi_plus[i] - qoi_minus[i]) / (2.*delta_p);
00767 
00768       this->assembly(true, false);
00769       this->rhs->close();
00770       *partialR_partialp += *this->rhs;
00771       *partialR_partialp /= (2.*delta_p);
00772 
00773       // Don't leave the parameter changed
00774       *parameters[j] = old_parameter;
00775 
00776       for (unsigned int i=0; i != Nq; ++i)
00777         if (qoi_indices.has_index(i))
00778           sensitivities[i][j] = partialq_partialp[i] -
00779             partialR_partialp->dot(this->get_adjoint_solution(i));
00780     }
00781 
00782   // All parameters have been reset.
00783   // We didn't cache the original rhs or matrix for memory reasons,
00784   // but we can restore them to a state consistent solution -
00785   // principle of least surprise.
00786   this->assembly(true, true);
00787   this->rhs->close();
00788   this->matrix->close();
00789   this->assemble_qoi(qoi_indices);
00790 }

std::pair< unsigned int, Real > libMesh::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 libMesh::System.

Reimplemented in libMesh::DifferentiableSystem.

Definition at line 382 of file implicit_system.C.

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

00383 {
00384   // Log how long the linear solve takes.
00385   START_LOG("adjoint_solve()", "ImplicitSystem");
00386 
00387   if (this->assemble_before_solve)
00388     // Assemble the linear system
00389     this->assembly (/* get_residual = */ false,
00390                     /* get_jacobian = */ true);
00391 
00392   // The adjoint problem is linear
00393   LinearSolver<Number> *linear_solver = this->get_linear_solver();
00394 
00395   // Reset and build the RHS from the QOI derivative
00396   this->assemble_qoi_derivative(qoi_indices);
00397 
00398   // Our iteration counts and residuals will be sums of the individual
00399   // results
00400   std::pair<unsigned int, Real> solver_params =
00401     this->get_linear_solve_parameters();
00402   std::pair<unsigned int, Real> totalrval = std::make_pair(0,0.0);
00403 
00404   for (unsigned int i=0; i != this->qoi.size(); ++i)
00405     if (qoi_indices.has_index(i))
00406       {
00407         const std::pair<unsigned int, Real> rval =
00408           linear_solver->adjoint_solve (*matrix, this->add_adjoint_solution(i),
00409                                         this->get_adjoint_rhs(i),
00410                                         solver_params.second,
00411                                         solver_params.first);
00412 
00413             totalrval.first  += rval.first;
00414             totalrval.second += rval.second;
00415       }
00416 
00417   this->release_linear_solver(linear_solver);
00418 
00419   // The linear solver may not have fit our constraints exactly
00420 #ifdef LIBMESH_ENABLE_CONSTRAINTS
00421   for (unsigned int i=0; i != this->qoi.size(); ++i)
00422     if (qoi_indices.has_index(i))
00423       this->get_dof_map().enforce_constraints_exactly
00424         (*this, &this->get_adjoint_solution(i),
00425          /* homogeneous = */ true);
00426 #endif
00427 
00428   // Stop logging the nonlinear solve
00429   STOP_LOG("adjoint_solve()", "ImplicitSystem");
00430 
00431   return totalrval;
00432 }

virtual void libMesh::LinearImplicitSystem::assemble ( ) [inline, virtual]

Prepares matrix and _dof_map for matrix assembly. Does not actually assemble anything. For matrix assembly, use the assemble() in derived classes. Should be overloaded in derived classes.

Reimplemented from libMesh::ImplicitSystem.

Reimplemented in libMesh::FrequencySystem, and libMesh::NewmarkSystem.

Definition at line 104 of file linear_implicit_system.h.

Referenced by assembly(), and solve().

00104 { ImplicitSystem::assemble(); }

void libMesh::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 libMesh::System.

Reimplemented in libMesh::FEMSystem.

Definition at line 89 of file explicit_system.C.

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

00090 {
00091   // The user quantity of interest assembly gets to expect to
00092   // accumulate on initially zero values
00093   for (unsigned int i=0; i != qoi.size(); ++i)
00094     if (qoi_indices.has_index(i))
00095       qoi[i] = 0;
00096 
00097   Parent::assemble_qoi (qoi_indices);
00098 }

void libMesh::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 libMesh::System.

Reimplemented in libMesh::FEMSystem.

Definition at line 102 of file explicit_system.C.

References libMesh::System::add_adjoint_rhs(), libMesh::System::assemble_qoi_derivative(), libMesh::QoISet::has_index(), libMesh::System::qoi, and libMesh::NumericVector< T >::zero().

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

00103 {
00104   // The user quantity of interest derivative assembly gets to expect
00105   // to accumulate on initially zero vectors
00106   for (unsigned int i=0; i != qoi.size(); ++i)
00107     if (qoi_indices.has_index(i))
00108       this->add_adjoint_rhs(i).zero();
00109 
00110   Parent::assemble_qoi_derivative (qoi_indices);
00111 }

void libMesh::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 libMesh::System.

Definition at line 660 of file implicit_system.C.

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

Referenced by libMesh::ImplicitSystem::sensitivity_solve().

00662 00663 00664 00665 00666 00667     { 00668 00669 00670 00671 00672 00673 00674 00675 00676 00677 00678 00679 00680 00681 00682 00683 00684 00685 00686 00687 00688 00689 00690     } 00691 }

class="fragment">00661 { const unsigned int Np = libmesh_cast_int<unsigned int> (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 libMesh::LinearImplicitSystem::assembly	(	bool	get_residual,
		bool	get_jacobian
	)			`[virtual]`

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

Reimplemented from libMesh::ImplicitSystem.

Definition at line 373 of file linear_implicit_system.C.

References libMesh::NumericVector< T >::add_vector(), assemble(), libMesh::SparseMatrix< T >::close(), libMesh::NumericVector< T >::close(), libMesh::ImplicitSystem::matrix, libMesh::ExplicitSystem::rhs, and libMesh::System::solution.

00375 {
00376   // Residual R(u(p),p) := A(p)*u(p) - b(p)
00377   // partial R / partial u = A
00378 
00379   this->assemble();
00380   this->rhs->close();
00381   this->matrix->close();
00382 
00383   *(this->rhs) *= -1.0;
00384   this->rhs->add_vector(*this->solution, *this->matrix);
00385 }

void libMesh::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 1753 of file system.C.

References libMesh::System::_assemble_system_function, libMesh::System::_assemble_system_object, and libMesh::out.

01755 {
01756   libmesh_assert(fptr);
01757 
01758   if (_assemble_system_object != NULL)
01759     {
01760       libmesh_here();
01761       libMesh::out << "WARNING:  Cannot specify both assembly function and object!"
01762                    << std::endl;
01763 
01764       _assemble_system_object = NULL;
01765     }
01766 
01767   _assemble_system_function = fptr;
01768 }

void libMesh::System::attach_assemble_object ( System::Assembly & assemble_in ) [inherited]

Register a user object to use in assembling the system matrix and RHS.

Definition at line 1772 of file system.C.

References libMesh::System::_assemble_system_function, libMesh::System::_assemble_system_object, and libMesh::out.

01773 {
01774   if (_assemble_system_function != NULL)
01775     {
01776       libmesh_here();
01777       libMesh::out << "WARNING:  Cannot specify both assembly object and function!"
01778                    << std::endl;
01779 
01780       _assemble_system_function = NULL;
01781     }
01782 
01783   _assemble_system_object = &assemble_in;
01784 }

void libMesh::System::attach_constraint_function ( void fptrEquationSystems &es,const std::string &name ) [inherited]

Register a user function for imposing constraints.

Definition at line 1788 of file system.C.

References libMesh::System::_constrain_system_function, libMesh::System::_constrain_system_object, and libMesh::out.

01790 {
01791   libmesh_assert(fptr);
01792 
01793   if (_constrain_system_object != NULL)
01794     {
01795       libmesh_here();
01796       libMesh::out << "WARNING:  Cannot specify both constraint function and object!"
01797                    << std::endl;
01798 
01799       _constrain_system_object = NULL;
01800     }
01801 
01802   _constrain_system_function = fptr;
01803 }

void libMesh::System::attach_constraint_object ( System::Constraint & constrain ) [inherited]

Register a user object for imposing constraints.

Definition at line 1807 of file system.C.

References libMesh::System::_constrain_system_function, libMesh::System::_constrain_system_object, and libMesh::out.

01808 {
01809   if (_constrain_system_function != NULL)
01810     {
01811       libmesh_here();
01812       libMesh::out << "WARNING:  Cannot specify both constraint object and function!"
01813                    << std::endl;
01814 
01815       _constrain_system_function = NULL;
01816     }
01817 
01818   _constrain_system_object = &constrain;
01819 }

void libMesh::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 1718 of file system.C.

References libMesh::System::_init_system_function, libMesh::System::_init_system_object, and libMesh::out.

01720 {
01721   libmesh_assert(fptr);
01722 
01723   if (_init_system_object != NULL)
01724     {
01725       libmesh_here();
01726       libMesh::out << "WARNING:  Cannot specify both initialization function and object!"
01727                    << std::endl;
01728 
01729       _init_system_object = NULL;
01730     }
01731 
01732   _init_system_function = fptr;
01733 }

void libMesh::System::attach_init_object ( System::Initialization & init_in ) [inherited]

Register a user class to use to initialize the system. Note this is exclusive with the attach_init_function.

Definition at line 1737 of file system.C.

References libMesh::System::_init_system_function, libMesh::System::_init_system_object, and libMesh::out.

01738 {
01739   if (_init_system_function != NULL)
01740     {
01741       libmesh_here();
01742       libMesh::out << "WARNING:  Cannot specify both initialization object and function!"
01743                    << std::endl;
01744 
01745       _init_system_function = NULL;
01746     }
01747 
01748   _init_system_object = &init_in;
01749 }

void libMesh::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 libMesh::System::attach_QOI_derivative_object ( QOIDerivative & qoi_derivative ) [inherited]

Register a user object for evaluating derivatives of a quantity of interest with respect to test functions, whose values should be placed in System::rhs

Definition at line 1879 of file system.C.

References libMesh::System::_qoi_evaluate_derivative_function, libMesh::System::_qoi_evaluate_derivative_object, and libMesh::out.

01880 {
01881   if (_qoi_evaluate_derivative_function != NULL)
01882     {
01883       libmesh_here();
01884       libMesh::out << "WARNING:  Cannot specify both QOI derivative object and function!"
01885                    << std::endl;
01886 
01887       _qoi_evaluate_derivative_function = NULL;
01888     }
01889 
01890   _qoi_evaluate_derivative_object = &qoi_derivative;
01891 }

void libMesh::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 libMesh::System::attach_QOI_object ( QOI & qoi ) [inherited]

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

Definition at line 1843 of file system.C.

References libMesh::System::_qoi_evaluate_function, libMesh::System::_qoi_evaluate_object, and libMesh::out.

01844 {
01845   if (_qoi_evaluate_function != NULL)
01846     {
01847       libmesh_here();
01848       libMesh::out << "WARNING:  Cannot specify both QOI object and function!"
01849                    << std::endl;
01850 
01851       _qoi_evaluate_function = NULL;
01852     }
01853 
01854   _qoi_evaluate_object = &qoi_in;
01855 }

void libMesh::LinearImplicitSystem::attach_shell_matrix ( ShellMatrix< Number > * shell_matrix )

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 165 of file linear_implicit_system.C.

References _shell_matrix.

Referenced by detach_shell_matrix().

00166 {
00167   _shell_matrix = shell_matrix;
00168 }

void libMesh::System::boundary_project_solution	(	const std::set< boundary_id_type > &	b,
		const std::vector< unsigned int > &	variables,
		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,
		const Parameters &	parameters
	)			`[inherited]`

Projects arbitrary boundary functions onto a vector of degree of freedom values for the current system. Only degrees of freedom which affect the function's trace on a boundary in the set b are affected. Only degrees of freedom associated with the variables listed in the vector variables are projected. The function value fptr and its gradient gptr are represented by function pointers. A gradient gptr is only required/used for projecting onto finite element spaces with continuous derivatives.

This method projects components of an arbitrary boundary function onto the solution via L2 projections and nodal interpolations on each element.

Definition at line 641 of file system_projection.C.

00652 {
00653   WrappedFunction<Number> f(*this, fptr, &parameters);
00654   WrappedFunction<Gradient> g(*this, gptr, &parameters);
00655   this->boundary_project_solution(b, variables, &f, &g);
00656 }

void libMesh::System::boundary_project_solution	(	const std::set< boundary_id_type > &	b,
		const std::vector< unsigned int > &	variables,
		FunctionBase< Number > *	f,
		FunctionBase< Gradient > *	g = `NULL`
	)			`[inherited]`

Projects arbitrary boundary functions onto a vector of degree of freedom values for the current system. Only degrees of freedom which affect the function's trace on a boundary in the set b are affected. Only degrees of freedom associated with the variables listed in the vector variables are projected. The function value f and its gradient g are user-provided cloneable functors. A gradient g is only required/used for projecting onto finite element spaces with continuous derivatives. If non-default Parameters are to be used, they can be provided in the parameters argument.

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

Definition at line 665 of file system_projection.C.

00669 {
00670   this->boundary_project_vector(b, variables, *solution, f, g);
00671 
00672   solution->localize(*current_local_solution);
00673 }

void libMesh::System::boundary_project_vector	(	const std::set< boundary_id_type > &	b,
		const std::vector< unsigned int > &	variables,
		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,
		const Parameters &	parameters,
		NumericVector< Number > &	new_vector
	)			const `[inherited]`

Projects arbitrary boundary functions onto a vector of degree of freedom values for the current system. Only degrees of freedom which affect the function's trace on a boundary in the set b are affected. Only degrees of freedom associated with the variables listed in the vector variables are projected. The function value fptr and its gradient gptr are represented by function pointers. A gradient gptr is only required/used for projecting onto finite element spaces with continuous derivatives.

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

Definition at line 684 of file system_projection.C.

00696 {
00697   WrappedFunction<Number> f(*this, fptr, &parameters);
00698   WrappedFunction<Gradient> g(*this, gptr, &parameters);
00699   this->boundary_project_vector(b, variables, new_vector, &f, &g);
00700 }

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

Projects arbitrary boundary functions onto a vector of degree of freedom values for the current system. Only degrees of freedom which affect the function's trace on a boundary in the set b are affected. Only degrees of freedom associated with the variables listed in the vector variables are projected. The function value f and its gradient g are user-provided cloneable functors. A gradient g is only required/used for projecting onto finite element spaces with continuous derivatives. If non-default Parameters are to be used, they can be provided in the parameters argument.

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

Definition at line 707 of file system_projection.C.

References libMesh::NumericVector< T >::close(), and libMesh::Threads::parallel_for().

00712 {
00713   START_LOG ("boundary_project_vector()", "System");
00714 
00715   Threads::parallel_for
00716     (ConstElemRange (this->get_mesh().active_local_elements_begin(),
00717                      this->get_mesh().active_local_elements_end() ),
00718      BoundaryProjectSolution(b, variables, *this, f, g,
00719                              this->get_equation_systems().parameters,
00720                              new_vector)
00721     );
00722 
00723   // We don't do SCALAR dofs when just projecting the boundary, so
00724   // we're done here.
00725 
00726   new_vector.close();
00727 
00728 #ifdef LIBMESH_ENABLE_CONSTRAINTS
00729   this->get_dof_map().enforce_constraints_exactly(*this, &new_vector);
00730 #endif
00731 
00732   STOP_LOG("boundary_project_vector()", "System");
00733 }

Real libMesh::System::calculate_norm	(	const 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 1400 of file system.C.

References libMesh::System::_dof_map, std::abs(), libMesh::MeshBase::active_local_elements_begin(), libMesh::MeshBase::active_local_elements_end(), libMesh::TypeTensor< T >::add_scaled(), libMesh::TypeVector< T >::add_scaled(), libMesh::FEGenericBase< Real >::build(), libMesh::CommWorld, libMesh::FEType::default_quadrature_rule(), libMeshEnums::DISCRETE_L1, libMeshEnums::DISCRETE_L2, libMeshEnums::DISCRETE_L_INF, libMesh::System::discrete_var_norm(), libMesh::DofMap::dof_indices(), libMesh::AutoPtr< Tp >::get(), libMesh::System::get_dof_map(), libMesh::System::get_mesh(), libMeshEnums::H1, libMeshEnums::H1_SEMINORM, libMeshEnums::H2, libMeshEnums::H2_SEMINORM, libMesh::SystemNorm::is_discrete(), libMeshEnums::L1, libMesh::NumericVector< T >::l1_norm(), libMeshEnums::L2, libMesh::NumericVector< T >::l2_norm(), libMeshEnums::L_INF, libMesh::NumericVector< T >::linfty_norm(), libMesh::NumericVector< T >::localize(), libMesh::Parallel::Communicator::max(), std::max(), libMesh::MeshBase::mesh_dimension(), libMesh::System::n_vars(), libMesh::TensorTools::norm_sq(), libMesh::Real, libMeshEnums::SERIAL, libMesh::TypeTensor< T >::size(), libMesh::TypeVector< T >::size(), libMesh::NumericVector< T >::size(), libMesh::TypeTensor< T >::size_sq(), libMesh::TypeVector< T >::size_sq(), libMesh::Parallel::Communicator::sum(), libMesh::SystemNorm::type(), libMesh::DofMap::variable_type(), libMeshEnums::W1_INF_SEMINORM, libMeshEnums::W2_INF_SEMINORM, libMesh::SystemNorm::weight(), and libMesh::SystemNorm::weight_sq().

01402 {
01403   // This function must be run on all processors at once
01404   parallel_only();
01405 
01406   START_LOG ("calculate_norm()", "System");
01407 
01408   // Zero the norm before summation
01409   Real v_norm = 0.;
01410 
01411   if (norm.is_discrete())
01412     {
01413       STOP_LOG ("calculate_norm()", "System");
01414       //Check to see if all weights are 1.0 and all types are equal
01415       FEMNormType norm_type0 = norm.type(0);
01416       unsigned int check_var = 0;
01417       for (; check_var != this->n_vars(); ++check_var)
01418         if((norm.weight(check_var) != 1.0) || (norm.type(check_var) != norm_type0))
01419           break;
01420 
01421       //All weights were 1.0 so just do the full vector discrete norm
01422       if(check_var == this->n_vars())
01423         {
01424           if(norm_type0 == DISCRETE_L1)
01425             return v.l1_norm();
01426           if(norm_type0 == DISCRETE_L2)
01427             return v.l2_norm();
01428           if(norm_type0 == DISCRETE_L_INF)
01429             return v.linfty_norm();
01430           else
01431             libmesh_error();
01432         }
01433 
01434       for (unsigned int var=0; var != this->n_vars(); ++var)
01435         {
01436           // Skip any variables we don't need to integrate
01437           if (norm.weight(var) == 0.0)
01438             continue;
01439 
01440           v_norm += norm.weight(var) * discrete_var_norm(v, var, norm.type(var));
01441         }
01442 
01443       return v_norm;
01444     }
01445 
01446   // Localize the potentially parallel vector
01447   AutoPtr<NumericVector<Number> > local_v = NumericVector<Number>::build();
01448   local_v->init(v.size(), true, SERIAL);
01449   v.localize (*local_v, _dof_map->get_send_list());
01450 
01451   unsigned int dim = this->get_mesh().mesh_dimension();
01452 
01453   // I'm not sure how best to mix Hilbert norms on some variables (for
01454   // which we'll want to square then sum then square root) with norms
01455   // like L_inf (for which we'll just want to take an absolute value
01456   // and then sum).
01457   bool using_hilbert_norm = true,
01458        using_nonhilbert_norm = true;
01459 
01460   // Loop over all variables
01461   for (unsigned int var=0; var != this->n_vars(); ++var)
01462     {
01463       // Skip any variables we don't need to integrate
01464       Real norm_weight_sq = norm.weight_sq(var);
01465       if (norm_weight_sq == 0.0)
01466         continue;
01467       Real norm_weight = norm.weight(var);
01468 
01469       // Check for unimplemented norms (rather than just returning 0).
01470       FEMNormType norm_type = norm.type(var);
01471       if((norm_type==H1) ||
01472          (norm_type==H2) ||
01473          (norm_type==L2) ||
01474          (norm_type==H1_SEMINORM) ||
01475          (norm_type==H2_SEMINORM))
01476         {
01477           if (!using_hilbert_norm)
01478             libmesh_not_implemented();
01479           using_nonhilbert_norm = false;
01480         }
01481       else if ((norm_type==L1) ||
01482                (norm_type==L_INF) ||
01483                (norm_type==W1_INF_SEMINORM) ||
01484                (norm_type==W2_INF_SEMINORM))
01485         {
01486           if (!using_nonhilbert_norm)
01487             libmesh_not_implemented();
01488           using_hilbert_norm = false;
01489         }
01490       else
01491         libmesh_not_implemented();
01492 
01493       const FEType& fe_type = this->get_dof_map().variable_type(var);
01494       AutoPtr<QBase> qrule =
01495         fe_type.default_quadrature_rule (dim);
01496       AutoPtr<FEBase> fe
01497         (FEBase::build(dim, fe_type));
01498       fe->attach_quadrature_rule (qrule.get());
01499 
01500       const std::vector<Real>&               JxW = fe->get_JxW();
01501       const std::vector<std::vector<Real> >* phi = NULL;
01502       if (norm_type == H1 ||
01503           norm_type == H2 ||
01504           norm_type == L2 ||
01505           norm_type == L1 ||
01506           norm_type == L_INF)
01507         phi = &(fe->get_phi());
01508 
01509       const std::vector<std::vector<RealGradient> >* dphi = NULL;
01510       if (norm_type == H1 ||
01511           norm_type == H2 ||
01512           norm_type == H1_SEMINORM ||
01513           norm_type == W1_INF_SEMINORM)
01514         dphi = &(fe->get_dphi());
01515 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
01516       const std::vector<std::vector<RealTensor> >*   d2phi = NULL;
01517       if (norm_type == H2 ||
01518           norm_type == H2_SEMINORM ||
01519           norm_type == W2_INF_SEMINORM)
01520         d2phi = &(fe->get_d2phi());
01521 #endif
01522 
01523       std::vector<dof_id_type> dof_indices;
01524 
01525       // Begin the loop over the elements
01526       MeshBase::const_element_iterator       el     =
01527         this->get_mesh().active_local_elements_begin();
01528       const MeshBase::const_element_iterator end_el =
01529         this->get_mesh().active_local_elements_end();
01530 
01531       for ( ; el != end_el; ++el)
01532         {
01533           const Elem* elem = *el;
01534 
01535           fe->reinit (elem);
01536 
01537           this->get_dof_map().dof_indices (elem, dof_indices, var);
01538 
01539           const unsigned int n_qp = qrule->n_points();
01540 
01541           const unsigned int n_sf = libmesh_cast_int<unsigned int>
01542             (dof_indices.size());
01543 
01544           // Begin the loop over the Quadrature points.
01545           for (unsigned int qp=0; qp<n_qp; qp++)
01546             {
01547               if (norm_type == L1)
01548                 {
01549                   Number u_h = 0.;
01550                   for (unsigned int i=0; i != n_sf; ++i)
01551                     u_h += (*phi)[i][qp] * (*local_v)(dof_indices[i]);
01552                   v_norm += norm_weight *
01553                             JxW[qp] * std::abs(u_h);
01554                 }
01555 
01556               if (norm_type == L_INF)
01557                 {
01558                   Number u_h = 0.;
01559                   for (unsigned int i=0; i != n_sf; ++i)
01560                     u_h += (*phi)[i][qp] * (*local_v)(dof_indices[i]);
01561                   v_norm = std::max(v_norm, norm_weight * std::abs(u_h));
01562                 }
01563 
01564               if (norm_type == H1 ||
01565                   norm_type == H2 ||
01566                   norm_type == L2)
01567                 {
01568                   Number u_h = 0.;
01569                   for (unsigned int i=0; i != n_sf; ++i)
01570                     u_h += (*phi)[i][qp] * (*local_v)(dof_indices[i]);
01571                   v_norm += norm_weight_sq *
01572                             JxW[qp] * TensorTools::norm_sq(u_h);
01573                 }
01574 
01575               if (norm_type == H1 ||
01576                   norm_type == H2 ||
01577                   norm_type == H1_SEMINORM)
01578                 {
01579                   Gradient grad_u_h;
01580                   for (unsigned int i=0; i != n_sf; ++i)
01581                     grad_u_h.add_scaled((*dphi)[i][qp], (*local_v)(dof_indices[i]));
01582                   v_norm += norm_weight_sq *
01583                             JxW[qp] * grad_u_h.size_sq();
01584                 }
01585 
01586               if (norm_type == W1_INF_SEMINORM)
01587                 {
01588                   Gradient grad_u_h;
01589                   for (unsigned int i=0; i != n_sf; ++i)
01590                     grad_u_h.add_scaled((*dphi)[i][qp], (*local_v)(dof_indices[i]));
01591                   v_norm = std::max(v_norm, norm_weight * grad_u_h.size());
01592                 }
01593 
01594 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
01595               if (norm_type == H2 ||
01596                   norm_type == H2_SEMINORM)
01597                 {
01598                   Tensor hess_u_h;
01599                   for (unsigned int i=0; i != n_sf; ++i)
01600                     hess_u_h.add_scaled((*d2phi)[i][qp], (*local_v)(dof_indices[i]));
01601                   v_norm += norm_weight_sq *
01602                             JxW[qp] * hess_u_h.size_sq();
01603                 }
01604 
01605               if (norm_type == W2_INF_SEMINORM)
01606                 {
01607                   Tensor hess_u_h;
01608                   for (unsigned int i=0; i != n_sf; ++i)
01609                     hess_u_h.add_scaled((*d2phi)[i][qp], (*local_v)(dof_indices[i]));
01610                   v_norm = std::max(v_norm, norm_weight * hess_u_h.size());
01611                 }
01612 #endif
01613             }
01614         }
01615     }
01616 
01617   if (using_hilbert_norm)
01618     {
01619       CommWorld.sum(v_norm);
01620       v_norm = std::sqrt(v_norm);
01621     }
01622   else
01623     {
01624       CommWorld.max(v_norm);
01625     }
01626 
01627   STOP_LOG ("calculate_norm()", "System");
01628 
01629   return v_norm;
01630 }

Real libMesh::System::calculate_norm	(	const 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, L_INF, H1)

Definition at line 1379 of file system.C.

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

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

01382 {
01383   //short circuit to save time
01384   if(norm_type == DISCRETE_L1 ||
01385      norm_type == DISCRETE_L2 ||
01386      norm_type == DISCRETE_L_INF)
01387     return discrete_var_norm(v,var,norm_type);
01388 
01389   // Not a discrete norm
01390   std::vector<FEMNormType> norms(this->n_vars(), L2);
01391   std::vector<Real> weights(this->n_vars(), 0.0);
01392   norms[var] = norm_type;
01393   weights[var] = 1.0;
01394   Real val = this->calculate_norm(v, SystemNorm(norms, weights));
01395   return val;
01396 }

void libMesh::LinearImplicitSystem::clear ( ) [virtual]

Clear all the data structures associated with the system.

Reimplemented from libMesh::ImplicitSystem.

Reimplemented in libMesh::FrequencySystem, libMesh::NewmarkSystem, libMesh::RBConstruction, libMesh::RBEIMConstruction, libMesh::TransientRBConstruction, libMesh::RBConstructionBase< LinearImplicitSystem >, and libMesh::TransientSystem< RBConstruction >.

Definition at line 61 of file linear_implicit_system.C.

References libMesh::ImplicitSystem::clear(), linear_solver, and restrict_solve_to().

Referenced by ~LinearImplicitSystem().

00062 {
00063   // clear the linear solver
00064   linear_solver->clear();
00065 
00066   this->restrict_solve_to(NULL);
00067 
00068   // clear the parent data
00069   Parent::clear();
00070 }

bool libMesh::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 525 of file system.C.

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

Referenced by libMesh::EquationSystems::compare().

00528 {
00529   // we do not care for matrices, but for vectors
00530   libmesh_assert (!_can_add_vectors);
00531   libmesh_assert (!other_system._can_add_vectors);
00532 
00533   if (verbose)
00534     {
00535       libMesh::out << "  Systems \"" << _sys_name << "\"" << std::endl;
00536       libMesh::out << "   comparing matrices not supported." << std::endl;
00537       libMesh::out << "   comparing names...";
00538     }
00539 
00540   // compare the name: 0 means identical
00541   const int name_result = _sys_name.compare(other_system.name());
00542   if (verbose)
00543     {
00544       if (name_result == 0)
00545         libMesh::out << " identical." << std::endl;
00546       else
00547         libMesh::out << "  names not identical." << std::endl;
00548       libMesh::out << "   comparing solution vector...";
00549     }
00550 
00551 
00552   // compare the solution: -1 means identical
00553   const int solu_result = solution->compare (*other_system.solution.get(),
00554                                              threshold);
00555 
00556   if (verbose)
00557     {
00558       if (solu_result == -1)
00559         libMesh::out << " identical up to threshold." << std::endl;
00560       else
00561         libMesh::out << "  first difference occured at index = "
00562                   << solu_result << "." << std::endl;
00563     }
00564 
00565 
00566   // safety check, whether we handle at least the same number
00567   // of vectors
00568   std::vector<int> ov_result;
00569 
00570   if (this->n_vectors() != other_system.n_vectors())
00571     {
00572       if (verbose)
00573         {
00574           libMesh::out << "   Fatal difference. This system handles "
00575                     << this->n_vectors() << " add'l vectors," << std::endl
00576                     << "   while the other system handles "
00577                     << other_system.n_vectors()
00578                     << " add'l vectors." << std::endl
00579                     << "   Aborting comparison." << std::endl;
00580         }
00581       return false;
00582     }
00583   else if (this->n_vectors() == 0)
00584     {
00585       // there are no additional vectors...
00586       ov_result.clear ();
00587     }
00588   else
00589     {
00590       // compare other vectors
00591       for (const_vectors_iterator pos = _vectors.begin();
00592            pos != _vectors.end(); ++pos)
00593         {
00594           if (verbose)
00595               libMesh::out << "   comparing vector \""
00596                         << pos->first << "\" ...";
00597 
00598           // assume they have the same name
00599           const NumericVector<Number>& other_system_vector =
00600               other_system.get_vector(pos->first);
00601 
00602           ov_result.push_back(pos->second->compare (other_system_vector,
00603                                                     threshold));
00604 
00605           if (verbose)
00606             {
00607               if (ov_result[ov_result.size()-1] == -1)
00608                 libMesh::out << " identical up to threshold." << std::endl;
00609               else
00610                 libMesh::out << " first difference occured at" << std::endl
00611                           << "   index = " << ov_result[ov_result.size()-1] << "." << std::endl;
00612             }
00613 
00614         }
00615 
00616     } // finished comparing additional vectors
00617 
00618 
00619   bool overall_result;
00620 
00621   // sum up the results
00622   if ((name_result==0) && (solu_result==-1))
00623     {
00624       if (ov_result.size()==0)
00625         overall_result = true;
00626       else
00627         {
00628           bool ov_identical;
00629           unsigned int n    = 0;
00630           do
00631             {
00632               ov_identical = (ov_result[n]==-1);
00633               n++;
00634             }
00635           while (ov_identical && n<ov_result.size());
00636           overall_result = ov_identical;
00637         }
00638     }
00639   else
00640     overall_result = false;
00641 
00642   if (verbose)
00643     {
00644       libMesh::out << "   finished comparisons, ";
00645       if (overall_result)
00646         libMesh::out << "found no differences." << std::endl << std::endl;
00647       else
00648         libMesh::out << "found differences." << std::endl << std::endl;
00649     }
00650 
00651   return overall_result;
00652 }

Number libMesh::System::current_solution ( const dof_id_type global_dof_number ) const [inherited]

Returns:: the current solution for the specified global DOF.

Definition at line 181 of file system.C.

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

Referenced by libMesh::ExactSolution::_compute_error(), libMesh::HPCoarsenTest::add_projection(), libMesh::JumpErrorEstimator::estimate_error(), libMesh::ExactErrorEstimator::estimate_error(), libMesh::WeightedPatchRecoveryErrorEstimator::EstimateError::operator()(), libMesh::PatchRecoveryErrorEstimator::EstimateError::operator()(), libMesh::System::point_gradient(), libMesh::System::point_hessian(), libMesh::System::point_value(), libMesh::FEMContext::pre_fe_reinit(), libMesh::HPCoarsenTest::select_refinement(), libMesh::EnsightIO::write_scalar_ascii(), and libMesh::EnsightIO::write_vector_ascii().

00182 {
00183   // Check the sizes
00184   libmesh_assert_less (global_dof_number, _dof_map->n_dofs());
00185   libmesh_assert_less (global_dof_number, current_local_solution->size());
00186 
00187   return (*current_local_solution)(global_dof_number);
00188 }

void libMesh::System::deactivate ( ) [inline, inherited]

Deactivates the system. Only active systems are solved.

Definition at line 1882 of file system.h.

References libMesh::System::_active.

01883 {
01884   _active = false;
01885 }

void libMesh::LinearImplicitSystem::detach_shell_matrix ( void ) [inline]

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

Definition at line 176 of file linear_implicit_system.h.

References attach_shell_matrix().

00176 { attach_shell_matrix(NULL); }

void libMesh::ImplicitSystem::disable_cache ( ) [virtual, inherited]

Assembles & solves the linear system Ax=b. Avoids use of any cached data that might affect any solve result. Should be overloaded in derived systems.

Reimplemented from libMesh::System.

Definition at line 313 of file implicit_system.C.

References libMesh::System::assemble_before_solve, libMesh::ImplicitSystem::get_linear_solver(), and libMesh::LinearSolver< T >::reuse_preconditioner().

00313                                     {
00314   this->assemble_before_solve = true;
00315   this->get_linear_solver()->reuse_preconditioner(false);
00316 }

void libMesh::ReferenceCounter::disable_print_counter_info ( ) [static, inherited]

Definition at line 106 of file reference_counter.C.

References libMesh::ReferenceCounter::_enable_print_counter.

00107 {
00108   _enable_print_counter = false;
00109   return;
00110 }

void libMesh::ReferenceCounter::enable_print_counter_info ( ) [static, inherited]

Methods to enable/disable the reference counter output from print_info()

Definition at line 100 of file reference_counter.C.

References libMesh::ReferenceCounter::_enable_print_counter.

00101 {
00102   _enable_print_counter = true;
00103   return;
00104 }

Real libMesh::LinearImplicitSystem::final_linear_residual ( ) const [inline]

Returns the final residual for the linear system solve.

Definition at line 160 of file linear_implicit_system.h.

References _final_linear_residual.

00160 { return _final_linear_residual; }

void libMesh::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 libMesh::System.

Definition at line 795 of file implicit_system.C.

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

00798 {
00799   const unsigned int Np = libmesh_cast_int<unsigned int>
00800     (parameters.size());
00801   const unsigned int Nq = libmesh_cast_int<unsigned int>
00802     (qoi.size());
00803 
00804   // We currently get partial derivatives via central differencing
00805   const Real delta_p = TOLERANCE;
00806 
00807   // An introduction to the problem:
00808   //
00809   // Residual R(u(p),p) = 0
00810   // partial R / partial u = J = system matrix
00811   //
00812   // This implies that:
00813   // d/dp(R) = 0
00814   // (partial R / partial p) +
00815   // (partial R / partial u) * (partial u / partial p) = 0
00816 
00817   // We first solve for (partial u / partial p) for each parameter:
00818   // J * (partial u / partial p) = - (partial R / partial p)
00819 
00820   this->sensitivity_solve(parameters);
00821 
00822   // Get ready to fill in senstivities:
00823   sensitivities.allocate_data(qoi_indices, *this, parameters);
00824 
00825   // We use the identity:
00826   // dq/dp = (partial q / partial p) + (partial q / partial u) *
00827   //         (partial u / partial p)
00828 
00829   // We get (partial q / partial u) from the user
00830   this->assemble_qoi_derivative(qoi_indices);
00831 
00832   // We don't need these to be closed() in this function, but libMesh
00833   // standard practice is to have them closed() by the time the
00834   // function exits
00835   for (unsigned int i=0; i != this->qoi.size(); ++i)
00836     if (qoi_indices.has_index(i))
00837       this->get_adjoint_rhs(i).close();
00838 
00839   for (unsigned int j=0; j != Np; ++j)
00840     {
00841       // (partial q / partial p) ~= (q(p+dp)-q(p-dp))/(2*dp)
00842 
00843       Number old_parameter = *parameters[j];
00844 
00845       *parameters[j] = old_parameter - delta_p;
00846       this->assemble_qoi();
00847       std::vector<Number> qoi_minus = this->qoi;
00848 
00849       *parameters[j] = old_parameter + delta_p;
00850       this->assemble_qoi();
00851       std::vector<Number>& qoi_plus = this->qoi;
00852 
00853       std::vector<Number> partialq_partialp(Nq, 0);
00854       for (unsigned int i=0; i != Nq; ++i)
00855         if (qoi_indices.has_index(i))
00856           partialq_partialp[i] = (qoi_plus[i] - qoi_minus[i]) / (2.*delta_p);
00857 
00858       // Don't leave the parameter changed
00859       *parameters[j] = old_parameter;
00860 
00861       for (unsigned int i=0; i != Nq; ++i)
00862         if (qoi_indices.has_index(i))
00863           sensitivities[i][j] = partialq_partialp[i] +
00864             this->get_adjoint_rhs(i).dot(this->get_sensitivity_solution(j));
00865     }
00866 
00867   // All parameters have been reset.
00868   // We didn't cache the original rhs or matrix for memory reasons,
00869   // but we can restore them to a state consistent solution -
00870   // principle of least surprise.
00871   this->assembly(true, true);
00872   this->rhs->close();
00873   this->matrix->close();
00874   this->assemble_qoi(qoi_indices);
00875 }

const NumericVector< Number > & libMesh::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 1036 of file system.C.

References libMesh::System::get_vector().

01037 {
01038   std::ostringstream adjoint_rhs_name;
01039   adjoint_rhs_name << "adjoint_rhs" << i;
01040 
01041   return this->get_vector(adjoint_rhs_name.str());
01042 }

NumericVector< Number > & libMesh::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 1026 of file system.C.

References libMesh::System::get_vector().

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

01027 {
01028   std::ostringstream adjoint_rhs_name;
01029   adjoint_rhs_name << "adjoint_rhs" << i;
01030 
01031   return this->get_vector(adjoint_rhs_name.str());
01032 }

const NumericVector< Number > & libMesh::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 976 of file system.C.

References libMesh::System::get_vector().

00977 {
00978   std::ostringstream adjoint_name;
00979   adjoint_name << "adjoint_solution" << i;
00980 
00981   return this->get_vector(adjoint_name.str());
00982 }

NumericVector< Number > & libMesh::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 966 of file system.C.

References libMesh::System::get_vector().

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

00967 {
00968   std::ostringstream adjoint_name;
00969   adjoint_name << "adjoint_solution" << i;
00970 
00971   return this->get_vector(adjoint_name.str());
00972 }

void libMesh::System::get_all_variable_numbers ( std::vector< unsigned int > & all_variable_numbers ) const [inherited]

Fills all_variable_numbers with all the variable numbers for the variables that have been added to this system.

Definition at line 1259 of file system.C.

References libMesh::System::_variable_numbers, and libMesh::System::n_vars().

01260 {
01261   all_variable_numbers.resize(n_vars());
01262 
01263   // Make sure the variable exists
01264   std::map<std::string, unsigned short int>::const_iterator
01265     it = _variable_numbers.begin();
01266   std::map<std::string, unsigned short int>::const_iterator
01267     it_end = _variable_numbers.end();
01268 
01269   unsigned int count = 0;
01270   for( ; it != it_end; ++it)
01271   {
01272     all_variable_numbers[count] = it->second;
01273     count++;
01274   }
01275 }

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

Returns:: a writeable reference to this system's _dof_map.

Definition at line 1858 of file system.h.

References libMesh::System::_dof_map.

01859 {
01860   return *_dof_map;
01861 }

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

Returns:: a constant reference to this system's _dof_map.

Definition at line 1850 of file system.h.

References libMesh::System::_dof_map.

Referenced by libMesh::__libmesh_petsc_diff_solver_jacobian(), libMesh::__libmesh_petsc_diff_solver_residual(), libMesh::ExactSolution::_compute_error(), libMesh::HPCoarsenTest::add_projection(), libMesh::UnsteadySolver::adjoint_advance_timestep(), libMesh::ImplicitSystem::adjoint_solve(), libMesh::UnsteadySolver::advance_timestep(), libMesh::EquationSystems::allgather(), libMesh::EquationSystems::build_discontinuous_solution_vector(), libMesh::EquationSystems::build_solution_vector(), libMesh::System::calculate_norm(), libMesh::Problem_Interface::computeF(), libMesh::Problem_Interface::computeJacobian(), libMesh::Problem_Interface::computePreconditioner(), DMCreateDomainDecomposition_libMesh(), DMCreateFieldDecomposition_libMesh(), DMFunction_libMesh(), DMJacobian_libMesh(), DMLibMeshSetSystem(), libMesh::DofMap::enforce_constraints_exactly(), libMesh::JumpErrorEstimator::estimate_error(), libMesh::ExactErrorEstimator::estimate_error(), libMesh::AdjointRefinementEstimator::estimate_error(), libMesh::System::get_info(), libMesh::EquationSystems::get_solution(), libMesh::SystemSubsetBySubdomain::init(), libMesh::UnsteadySolver::init_data(), libMesh::ImplicitSystem::init_matrices(), libMesh::EigenSystem::init_matrices(), libMesh::CondensedEigenSystem::initialize_condensed_dofs(), libMesh::System::local_dof_indices(), libMesh::DofMap::max_constraint_error(), libMesh::UnsteadySolver::old_nonlinear_solution(), libMesh::WeightedPatchRecoveryErrorEstimator::EstimateError::operator()(), libMesh::BoundaryProjectSolution::operator()(), libMesh::ProjectFEMSolution::operator()(), libMesh::ProjectSolution::operator()(), libMesh::PatchRecoveryErrorEstimator::EstimateError::operator()(), libMesh::ErrorVector::plot_error(), libMesh::System::point_gradient(), libMesh::System::point_hessian(), libMesh::System::point_value(), libMesh::FEMContext::pre_fe_reinit(), libMesh::System::project_vector(), libMesh::System::re_update(), libMesh::System::read_parallel_data(), libMesh::System::read_SCALAR_dofs(), libMesh::UnsteadySolver::reinit(), libMesh::ImplicitSystem::reinit(), libMesh::EquationSystems::reinit(), libMesh::EigenSystem::reinit(), libMesh::UnsteadySolver::retrieve_timestep(), libMesh::HPCoarsenTest::select_refinement(), libMesh::ImplicitSystem::sensitivity_solve(), libMesh::PetscDiffSolver::solve(), libMesh::NewtonSolver::solve(), libMesh::ImplicitSystem::weighted_sensitivity_adjoint_solve(), libMesh::ImplicitSystem::weighted_sensitivity_solve(), libMesh::System::write_parallel_data(), libMesh::EnsightIO::write_scalar_ascii(), libMesh::System::write_SCALAR_dofs(), and libMesh::EnsightIO::write_vector_ascii().

01851 {
01852   return *_dof_map;
01853 }

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

Returns:: a reference to this system's parent EquationSystems object.

Definition at line 684 of file system.h.

References libMesh::System::_equation_systems.

00684 { return _equation_systems; }

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

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

Definition at line 679 of file system.h.

References libMesh::System::_equation_systems.

00679 { return _equation_systems; }

std::string libMesh::ReferenceCounter::get_info ( ) [static, inherited]

Gets a string containing the reference information.

Definition at line 47 of file reference_counter.C.

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

Referenced by libMesh::ReferenceCounter::print_info().

00048 {
00049 #if defined(LIBMESH_ENABLE_REFERENCE_COUNTING) && defined(DEBUG)
00050 
00051   std::ostringstream oss;
00052 
00053   oss << '\n'
00054       << " ---------------------------------------------------------------------------- \n"
00055       << "| Reference count information                                                |\n"
00056       << " ---------------------------------------------------------------------------- \n";
00057 
00058   for (Counts::iterator it = _counts.begin();
00059        it != _counts.end(); ++it)
00060     {
00061       const std::string name(it->first);
00062       const unsigned int creations    = it->second.first;
00063       const unsigned int destructions = it->second.second;
00064 
00065       oss << "| " << name << " reference count information:\n"
00066           << "|  Creations:    " << creations    << '\n'
00067           << "|  Destructions: " << destructions << '\n';
00068     }
00069 
00070   oss << " ---------------------------------------------------------------------------- \n";
00071 
00072   return oss.str();
00073 
00074 #else
00075 
00076   return "";
00077 
00078 #endif
00079 }

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

Returns:: a string containing information about the system.

Definition at line 1634 of file system.C.

References libMesh::Utility::enum_to_string< FEFamily >(), libMesh::Utility::enum_to_string< InfMapType >(), libMesh::Utility::enum_to_string< Order >(), libMesh::FEType::family, libMesh::System::get_dof_map(), libMesh::DofMap::get_info(), libMesh::FEType::inf_map, libMesh::System::n_constrained_dofs(), libMesh::System::n_dofs(), libMesh::System::n_local_constrained_dofs(), libMesh::System::n_local_dofs(), libMesh::System::n_matrices(), libMesh::System::n_variable_groups(), libMesh::VariableGroup::n_variables(), libMesh::System::n_vectors(), libMesh::VariableGroup::name(), libMesh::System::name(), libMesh::System::number(), libMesh::FEType::order, libMesh::FEType::radial_family, libMesh::FEType::radial_order, libMesh::System::system_type(), libMesh::Variable::type(), libMesh::DofMap::variable_group(), and libMesh::System::variable_group().

01635 {
01636   std::ostringstream oss;
01637 
01638 
01639   const std::string& sys_name = this->name();
01640 
01641   oss << "   System #"  << this->number() << ", \"" << sys_name << "\"\n"
01642       << "    Type \""  << this->system_type() << "\"\n"
01643       << "    Variables=";
01644 
01645   for (unsigned int vg=0; vg<this->n_variable_groups(); vg++)
01646     {
01647       const VariableGroup &vg_description (this->variable_group(vg));
01648 
01649       if (vg_description.n_variables() > 1) oss << "{ ";
01650       for (unsigned int vn=0; vn<vg_description.n_variables(); vn++)
01651         oss << "\"" << vg_description.name(vn) << "\" ";
01652       if (vg_description.n_variables() > 1) oss << "} ";
01653     }
01654 
01655   oss << '\n';
01656 
01657   oss << "    Finite Element Types=";
01658 #ifndef LIBMESH_ENABLE_INFINITE_ELEMENTS
01659   for (unsigned int vg=0; vg<this->n_variable_groups(); vg++)
01660     oss << "\""
01661         << Utility::enum_to_string<FEFamily>(this->get_dof_map().variable_group(vg).type().family)
01662         << "\" ";
01663 #else
01664   for (unsigned int vg=0; vg<this->n_variable_groups(); vg++)
01665     {
01666       oss << "\""
01667           << Utility::enum_to_string<FEFamily>(this->get_dof_map().variable_group(vg).type().family)
01668           << "\", \""
01669           << Utility::enum_to_string<FEFamily>(this->get_dof_map().variable_group(vg).type().radial_family)
01670           << "\" ";
01671     }
01672 
01673   oss << '\n' << "    Infinite Element Mapping=";
01674   for (unsigned int vg=0; vg<this->n_variable_groups(); vg++)
01675     oss << "\""
01676         << Utility::enum_to_string<InfMapType>(this->get_dof_map().variable_group(vg).type().inf_map)
01677         << "\" ";
01678 #endif
01679 
01680   oss << '\n';
01681 
01682   oss << "    Approximation Orders=";
01683   for (unsigned int vg=0; vg<this->n_variable_groups(); vg++)
01684     {
01685 #ifndef LIBMESH_ENABLE_INFINITE_ELEMENTS
01686       oss << "\""
01687           << Utility::enum_to_string<Order>(this->get_dof_map().variable_group(vg).type().order)
01688           << "\" ";
01689 #else
01690       oss << "\""
01691           << Utility::enum_to_string<Order>(this->get_dof_map().variable_group(vg).type().order)
01692           << "\", \""
01693           << Utility::enum_to_string<Order>(this->get_dof_map().variable_group(vg).type().radial_order)
01694           << "\" ";
01695 #endif
01696     }
01697 
01698   oss << '\n';
01699 
01700   oss << "    n_dofs()="             << this->n_dofs()             << '\n';
01701   oss << "    n_local_dofs()="       << this->n_local_dofs()       << '\n';
01702 #ifdef LIBMESH_ENABLE_CONSTRAINTS
01703   oss << "    n_constrained_dofs()=" << this->n_constrained_dofs() << '\n';
01704   oss << "    n_local_constrained_dofs()=" << this->n_local_constrained_dofs() << '\n';
01705 #endif
01706 
01707   oss << "    " << "n_vectors()="  << this->n_vectors()  << '\n';
01708   oss << "    " << "n_matrices()="  << this->n_matrices()  << '\n';
01709 //   oss << "    " << "n_additional_matrices()=" << this->n_additional_matrices() << '\n';
01710 
01711   oss << this->get_dof_map().get_info();
01712 
01713   return oss.str();
01714 }

std::pair< unsigned int, Real > libMesh::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 libMesh::DifferentiableSystem, and libMesh::NonlinearImplicitSystem.

Definition at line 1359 of file implicit_system.C.

References libMesh::System::get_equation_systems(), and libMesh::Real.

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

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

LinearSolver< Number > * libMesh::LinearImplicitSystem::get_linear_solver ( ) const [virtual]

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

Reimplemented from libMesh::ImplicitSystem.

Definition at line 360 of file linear_implicit_system.C.

References libMesh::AutoPtr< Tp >::get(), and linear_solver.

00361 {
00362   return linear_solver.get();
00363 }

SparseMatrix< Number > & libMesh::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 277 of file implicit_system.C.

References libMesh::ImplicitSystem::_matrices, and libMesh::err.

00278 {
00279   // Make sure the matrix exists
00280   matrices_iterator pos = _matrices.find (mat_name);
00281 
00282   if (pos == _matrices.end())
00283     {
00284       libMesh::err << "ERROR: matrix "
00285                     << mat_name
00286                     << " does not exist in this system!"
00287                     << std::endl;
00288       libmesh_error();
00289     }
00290 
00291   return *(pos->second);
00292 }

const SparseMatrix< Number > & libMesh::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 258 of file implicit_system.C.

References libMesh::ImplicitSystem::_matrices, and libMesh::err.

Referenced by libMesh::NewmarkSystem::compute_matrix(), libMesh::EigenTimeSolver::solve(), and libMesh::NewmarkSystem::update_rhs().

00259 {
00260   // Make sure the matrix exists
00261   const_matrices_iterator pos = _matrices.find (mat_name);
00262 
00263   if (pos == _matrices.end())
00264     {
00265       libMesh::err << "ERROR: matrix "
00266                     << mat_name
00267                     << " does not exist in this system!"
00268                     << std::endl;
00269       libmesh_error();
00270     }
00271 
00272   return *(pos->second);
00273 }

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

Returns:: a reference to this systems's _mesh.

Definition at line 1842 of file system.h.

References libMesh::System::_mesh.

01843 {
01844   return _mesh;
01845 }

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

Returns:: a constant reference to this systems's _mesh.

Definition at line 1834 of file system.h.

References libMesh::System::_mesh.

01835 {
01836   return _mesh;
01837 }

const NumericVector< Number > & libMesh::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 1066 of file system.C.

References libMesh::System::get_vector().

01067 {
01068   std::ostringstream sensitivity_rhs_name;
01069   sensitivity_rhs_name << "sensitivity_rhs" << i;
01070 
01071   return this->get_vector(sensitivity_rhs_name.str());
01072 }

NumericVector< Number > & libMesh::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 1056 of file system.C.

References libMesh::System::get_vector().

Referenced by libMesh::ImplicitSystem::sensitivity_solve().

01057 {
01058   std::ostringstream sensitivity_rhs_name;
01059   sensitivity_rhs_name << "sensitivity_rhs" << i;
01060 
01061   return this->get_vector(sensitivity_rhs_name.str());
01062 }

const NumericVector< Number > & libMesh::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 925 of file system.C.

References libMesh::System::get_vector().

00926 {
00927   std::ostringstream sensitivity_name;
00928   sensitivity_name << "sensitivity_solution" << i;
00929 
00930   return this->get_vector(sensitivity_name.str());
00931 }

NumericVector< Number > & libMesh::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 915 of file system.C.

References libMesh::System::get_vector().

Referenced by libMesh::ImplicitSystem::sensitivity_solve().

00916 {
00917   std::ostringstream sensitivity_name;
00918   sensitivity_name << "sensitivity_solution" << i;
00919 
00920   return this->get_vector(sensitivity_name.str());
00921 }

ShellMatrix<Number>* libMesh::LinearImplicitSystem::get_shell_matrix ( void ) [inline]

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

Definition at line 182 of file linear_implicit_system.h.

References _shell_matrix.

00182 { return _shell_matrix; }

NumericVector< Number > & libMesh::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 838 of file system.C.

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

00839 {
00840   vectors_iterator v = vectors_begin();
00841   vectors_iterator v_end = vectors_end();
00842   unsigned int num = 0;
00843   while((num<vec_num) && (v!=v_end))
00844     {
00845       num++;
00846       ++v;
00847     }
00848   libmesh_assert (v != v_end);
00849   return *(v->second);
00850 }

const NumericVector< Number > & libMesh::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 822 of file system.C.

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

00823 {
00824   const_vectors_iterator v = vectors_begin();
00825   const_vectors_iterator v_end = vectors_end();
00826   unsigned int num = 0;
00827   while((num<vec_num) && (v!=v_end))
00828     {
00829       num++;
00830       ++v;
00831     }
00832   libmesh_assert (v != v_end);
00833   return *(v->second);
00834 }

NumericVector< Number > & libMesh::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 803 of file system.C.

References libMesh::System::_vectors, and libMesh::err.

00804 {
00805   // Make sure the vector exists
00806   vectors_iterator pos = _vectors.find(vec_name);
00807 
00808   if (pos == _vectors.end())
00809     {
00810       libMesh::err << "ERROR: vector "
00811                     << vec_name
00812                     << " does not exist in this system!"
00813                     << std::endl;
00814       libmesh_error();
00815     }
00816 
00817   return *(pos->second);
00818 }

const NumericVector< Number > & libMesh::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 784 of file system.C.

References libMesh::System::_vectors, and libMesh::err.

Referenced by libMesh::UnsteadySolver::adjoint_advance_timestep(), libMesh::UnsteadySolver::advance_timestep(), libMesh::AdaptiveTimeSolver::advance_timestep(), libMesh::System::compare(), libMesh::UnsteadySolver::du(), libMesh::AdjointRefinementEstimator::estimate_error(), libMesh::System::get_adjoint_rhs(), libMesh::System::get_adjoint_solution(), libMesh::System::get_sensitivity_rhs(), libMesh::System::get_sensitivity_solution(), libMesh::System::get_weighted_sensitivity_adjoint_solution(), libMesh::System::get_weighted_sensitivity_solution(), libMesh::NewmarkSystem::initial_conditions(), libMesh::MemorySolutionHistory::retrieve(), libMesh::UnsteadySolver::retrieve_timestep(), libMesh::TwostepTimeSolver::solve(), libMesh::FrequencySystem::solve(), libMesh::NewmarkSystem::update_rhs(), and libMesh::NewmarkSystem::update_u_v_a().

00785 {
00786   // Make sure the vector exists
00787   const_vectors_iterator pos = _vectors.find(vec_name);
00788 
00789   if (pos == _vectors.end())
00790     {
00791       libMesh::err << "ERROR: vector "
00792                     << vec_name
00793                     << " does not exist in this system!"
00794                     << std::endl;
00795       libmesh_error();
00796     }
00797 
00798   return *(pos->second);
00799 }

const NumericVector< Number > & libMesh::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 1006 of file system.C.

References libMesh::System::get_vector().

01007 {
01008   std::ostringstream adjoint_name;
01009   adjoint_name << "weighted_sensitivity_adjoint_solution" << i;
01010 
01011   return this->get_vector(adjoint_name.str());
01012 }

NumericVector< Number > & libMesh::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 996 of file system.C.

References libMesh::System::get_vector().

Referenced by libMesh::ImplicitSystem::weighted_sensitivity_adjoint_solve().

00997 {
00998   std::ostringstream adjoint_name;
00999   adjoint_name << "weighted_sensitivity_adjoint_solution" << i;
01000 
01001   return this->get_vector(adjoint_name.str());
01002 }

const NumericVector< Number > & libMesh::System::get_weighted_sensitivity_solution ( ) const [inherited]

Returns:: a reference to the solution of the last weighted sensitivity solve

Definition at line 949 of file system.C.

References libMesh::System::get_vector().

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

NumericVector< Number > & libMesh::System::get_weighted_sensitivity_solution ( ) [inherited]

Returns:: a reference to the solution of the last weighted sensitivity solve

Definition at line 942 of file system.C.

References libMesh::System::get_vector().

Referenced by libMesh::ImplicitSystem::weighted_sensitivity_solve().

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

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

Returns:: true if a variable named var exists in this System

Definition at line 1232 of file system.C.

References libMesh::System::_variable_numbers.

Referenced by libMesh::GMVIO::copy_nodal_solution().

01233 {
01234   return _variable_numbers.count(var);
01235 }

bool libMesh::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 378 of file implicit_system.h.

References libMesh::ImplicitSystem::_matrices.

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

00379 {
00380   return (_matrices.count(mat_name));
00381 }

bool libMesh::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 2018 of file system.h.

References libMesh::System::_vectors.

Referenced by libMesh::System::add_vector(), and libMesh::System::remove_vector().

02019 {
02020   return (_vectors.count(vec_name));
02021 }

void libMesh::System::identify_variable_groups ( const bool ivg ) [inline, inherited]

Toggle automatic VariableGroup identification.

Definition at line 2002 of file system.h.

References libMesh::System::_identify_variable_groups.

02003 {
02004   _identify_variable_groups = ivg;
02005 }

bool libMesh::System::identify_variable_groups ( ) const [inline, inherited]

Returns:: true when VariableGroup structures should be automatically identified, false otherwise.

Definition at line 1994 of file system.h.

References libMesh::System::_identify_variable_groups.

Referenced by libMesh::System::add_variable().

01995 {
01996   return _identify_variable_groups;
01997 }

void libMesh::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 163 of file reference_counter.h.

References libMesh::ReferenceCounter::_counts, and libMesh::Threads::spin_mtx.

Referenced by libMesh::ReferenceCountedObject< RBParametrized >::ReferenceCountedObject().

00164 {
00165   Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
00166   std::pair<unsigned int, unsigned int>& p = _counts[name];
00167 
00168   p.first++;
00169 }

void libMesh::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 176 of file reference_counter.h.

References libMesh::ReferenceCounter::_counts, and libMesh::Threads::spin_mtx.

Referenced by libMesh::ReferenceCountedObject< RBParametrized >::~ReferenceCountedObject().

00177 {
00178   Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
00179   std::pair<unsigned int, unsigned int>& p = _counts[name];
00180 
00181   p.second++;
00182 }

void libMesh::System::init ( ) [inherited]

Initializes degrees of freedom on the current mesh. Sets the

Definition at line 223 of file system.C.

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

00224 {
00225   // First initialize any required data:
00226   // either only the basic System data
00227   if (_basic_system_only)
00228     System::init_data();
00229   // or all the derived class' data too
00230   else
00231     this->init_data();
00232 
00233   // If no variables have been added to this system
00234   // don't do anything
00235   if(!this->n_vars())
00236     return;
00237 
00238   // Then call the user-provided intialization function
00239   this->user_initialization();
00240 }

void libMesh::LinearImplicitSystem::init_data ( ) [virtual]

Initializes new data members of the system

Reimplemented from libMesh::ImplicitSystem.

Reimplemented in libMesh::FrequencySystem, libMesh::RBEIMConstruction, libMesh::RBConstructionBase< LinearImplicitSystem >, and libMesh::TransientSystem< RBConstruction >.

Definition at line 74 of file linear_implicit_system.C.

References libMesh::ImplicitSystem::init_data(), and linear_solver.

00075 {
00076   // initialize parent data
00077   Parent::init_data();
00078 
00079   // re-initialize the linear solver interface
00080   linear_solver->clear();
00081 }

void libMesh::ImplicitSystem::init_matrices ( ) [protected, virtual, inherited]

Initializes the matrices associated with this system.

Definition at line 105 of file implicit_system.C.

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

Referenced by libMesh::ImplicitSystem::init_data().

00106 {
00107   libmesh_assert(matrix);
00108 
00109   // Check for quick return in case the system matrix
00110   // (and by extension all the matrices) has already
00111   // been initialized
00112   if (matrix->initialized())
00113     return;
00114 
00115   // Get a reference to the DofMap
00116   DofMap& dof_map = this->get_dof_map();
00117 
00118   // no chance to add other matrices
00119   _can_add_matrices = false;
00120 
00121   // Tell the matrices about the dof map, and vice versa
00122   for (matrices_iterator pos = _matrices.begin();
00123        pos != _matrices.end(); ++pos)
00124     {
00125       SparseMatrix<Number> &m = *(pos->second);
00126       libmesh_assert (!m.initialized());
00127 
00128       // We want to allow repeated init() on systems, but we don't
00129       // want to attach the same matrix to the DofMap twice
00130       if (!dof_map.is_attached(m))
00131         dof_map.attach_matrix (m);
00132     }
00133 
00134   // Compute the sparsity pattern for the current
00135   // mesh and DOF distribution.  This also updates
00136   // additional matrices, \p DofMap now knows them
00137   dof_map.compute_sparsity (this->get_mesh());
00138 
00139   // Initialize matrices
00140   for (matrices_iterator pos = _matrices.begin();
00141        pos != _matrices.end(); ++pos)
00142     pos->second->init ();
00143 
00144   // Set the additional matrices to 0.
00145   for (matrices_iterator pos = _matrices.begin();
00146        pos != _matrices.end(); ++pos)
00147     pos->second->zero ();
00148 }

bool libMesh::System::is_adjoint_already_solved ( ) const [inline, inherited]

Accessor for the adjoint_already_solved boolean

Definition at line 359 of file system.h.

References libMesh::System::adjoint_already_solved.

Referenced by libMesh::AdjointResidualErrorEstimator::estimate_error().

00360   { return adjoint_already_solved;}

void libMesh::System::local_dof_indices	(	const unsigned int	var,
		std::set< dof_id_type > &	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 1278 of file system.C.

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

Referenced by libMesh::System::discrete_var_norm().

01280 {
01281   // Make sure the set is clear
01282   var_indices.clear();
01283 
01284   std::vector<dof_id_type> dof_indices;
01285 
01286   // Begin the loop over the elements
01287   MeshBase::const_element_iterator       el     =
01288     this->get_mesh().active_local_elements_begin();
01289   const MeshBase::const_element_iterator end_el =
01290     this->get_mesh().active_local_elements_end();
01291 
01292   const dof_id_type
01293     first_local = this->get_dof_map().first_dof(),
01294     end_local   = this->get_dof_map().end_dof();
01295 
01296   for ( ; el != end_el; ++el)
01297     {
01298       const Elem* elem = *el;
01299       this->get_dof_map().dof_indices (elem, dof_indices, var);
01300 
01301       for(unsigned int i=0; i<dof_indices.size(); i++)
01302         {
01303           dof_id_type dof = dof_indices[i];
01304 
01305           //If the dof is owned by the local processor
01306           if(first_local <= dof && dof < end_local)
01307             var_indices.insert(dof_indices[i]);
01308         }
01309     }
01310 }

dof_id_type libMesh::System::n_active_dofs ( ) const [inline, inherited]

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

Definition at line 2010 of file system.h.

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

02011 {
02012   return this->n_dofs() - this->n_constrained_dofs();
02013 }

unsigned int libMesh::System::n_components ( ) const [inline, inherited]

Returns:: the total number of scalar components in the system's variables. This will equal n_vars() in the case of all scalar-valued variables.

Definition at line 1914 of file system.h.

References libMesh::System::_variables, libMesh::Variable::first_scalar_number(), and libMesh::Variable::n_components().

Referenced by libMesh::System::add_variables(), libMesh::WrappedFunction< Output >::operator()(), and libMesh::System::project_vector().

01915 {
01916   if (_variables.empty())
01917     return 0;
01918 
01919   const Variable& last = _variables.back();
01920   return last.first_scalar_number() + last.n_components();
01921 }

dof_id_type libMesh::System::n_constrained_dofs ( ) const [inherited]

Returns:: the total number of constrained degrees of freedom in the system.

Definition at line 144 of file system.C.

References libMesh::System::_dof_map.

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

00145 {
00146 #ifdef LIBMESH_ENABLE_CONSTRAINTS
00147 
00148   return _dof_map->n_constrained_dofs();
00149 
00150 #else
00151 
00152   return 0;
00153 
00154 #endif
00155 }

dof_id_type libMesh::System::n_dofs ( ) const [inherited]

Returns:: the number of degrees of freedom in the system

Definition at line 137 of file system.C.

References libMesh::System::_dof_map.

Referenced by libMesh::System::add_vector(), libMesh::AdjointRefinementEstimator::estimate_error(), libMesh::System::get_info(), libMesh::UnsteadySolver::init_data(), libMesh::System::init_data(), libMesh::System::n_active_dofs(), libMesh::CondensedEigenSystem::n_global_non_condensed_dofs(), libMesh::System::project_vector(), libMesh::System::read_legacy_data(), libMesh::UnsteadySolver::reinit(), libMesh::System::reinit(), and libMesh::System::restrict_vectors().

00138 {
00139   return _dof_map->n_dofs();
00140 }

unsigned int libMesh::LinearImplicitSystem::n_linear_iterations ( ) const [inline]

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

Definition at line 155 of file linear_implicit_system.h.

References _n_linear_iterations.

00155 { return _n_linear_iterations; }

dof_id_type libMesh::System::n_local_constrained_dofs ( ) const [inherited]

Returns:: the number of constrained degrees of freedom on this processor.

Definition at line 159 of file system.C.

References libMesh::System::_dof_map.

Referenced by libMesh::System::get_info().

00160 {
00161 #ifdef LIBMESH_ENABLE_CONSTRAINTS
00162 
00163   return _dof_map->n_local_constrained_dofs();
00164 
00165 #else
00166 
00167   return 0;
00168 
00169 #endif
00170 }

dof_id_type libMesh::System::n_local_dofs ( ) const [inherited]

Returns:: the number of degrees of freedom local to this processor

Definition at line 174 of file system.C.

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

Referenced by libMesh::System::add_vector(), libMesh::AdjointRefinementEstimator::estimate_error(), libMesh::System::get_info(), libMesh::UnsteadySolver::init_data(), libMesh::System::init_data(), libMesh::System::project_vector(), libMesh::UnsteadySolver::reinit(), libMesh::System::reinit(), and libMesh::System::restrict_vectors().

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

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

Returns:: the number of matrices handled by this system

Reimplemented from libMesh::System.

Definition at line 385 of file implicit_system.h.

References libMesh::ImplicitSystem::_matrices.

00386 {
00387  return libmesh_cast_int<unsigned int>(_matrices.size());
00388 }

static unsigned int libMesh::ReferenceCounter::n_objects ( ) [inline, static, inherited]

Prints the number of outstanding (created, but not yet destroyed) objects.

Definition at line 79 of file reference_counter.h.

References libMesh::ReferenceCounter::_n_objects.

00080   { return _n_objects; }

unsigned int libMesh::System::n_variable_groups ( ) const [inline, inherited]

Returns:: the number of VariableGroup variable groups in the system

Definition at line 1906 of file system.h.

References libMesh::System::_variable_groups.

Referenced by libMesh::System::add_variable(), libMesh::System::get_info(), and libMesh::System::init_data().

01907 {
01908   return libmesh_cast_int<unsigned int>(_variable_groups.size());
01909 }

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

Returns:: the number of variables in the system

Definition at line 1898 of file system.h.

References libMesh::System::_variables.

Referenced by libMesh::DiffContext::add_localized_vector(), libMesh::System::add_variable(), libMesh::System::add_variables(), libMesh::EquationSystems::build_discontinuous_solution_vector(), libMesh::EquationSystems::build_solution_vector(), libMesh::System::calculate_norm(), libMesh::WrappedFunction< Output >::component(), libMesh::DiffContext::DiffContext(), libMesh::JumpErrorEstimator::estimate_error(), libMesh::ExactErrorEstimator::estimate_error(), libMesh::AdjointResidualErrorEstimator::estimate_error(), libMesh::ErrorEstimator::estimate_errors(), libMesh::ExactSolution::ExactSolution(), libMesh::FEMContext::FEMContext(), libMesh::System::get_all_variable_numbers(), libMesh::EquationSystems::get_solution(), libMesh::System::init(), libMesh::FEMSystem::init_context(), libMesh::WrappedFunction< Output >::operator()(), libMesh::WeightedPatchRecoveryErrorEstimator::EstimateError::operator()(), libMesh::ProjectSolution::operator()(), libMesh::PatchRecoveryErrorEstimator::EstimateError::operator()(), libMesh::FEMContext::pre_fe_reinit(), libMesh::System::project_vector(), libMesh::System::re_update(), libMesh::System::read_legacy_data(), libMesh::System::read_parallel_data(), libMesh::System::read_serialized_blocked_dof_objects(), libMesh::System::read_serialized_vector(), libMesh::System::read_serialized_vectors(), libMesh::System::reinit(), libMesh::EquationSystems::reinit(), libMesh::HPCoarsenTest::select_refinement(), libMesh::SystemSubsetBySubdomain::set_var_nums(), libMesh::System::write_header(), libMesh::System::write_parallel_data(), libMesh::System::write_serialized_blocked_dof_objects(), libMesh::System::write_serialized_vector(), libMesh::System::write_serialized_vectors(), and libMesh::System::zero_variable().

01899 {
01900   return libmesh_cast_int<unsigned int>(_variables.size());
01901 }

unsigned int libMesh::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 2026 of file system.h.

References libMesh::System::_vectors.

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

02027 {
02028   return libmesh_cast_int<unsigned int>(_vectors.size());
02029 }

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

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

Returns:: the system number.

Definition at line 1826 of file system.h.

References libMesh::System::_sys_number.

01827 {
01828   return _sys_number;
01829 }

Gradient libMesh::System::point_gradient	(	unsigned int	var,
		const Point &	p,
		const Elem &	e
	)			const `[inherited]`

Returns the gradient of the solution variable var at the physical point p in local Elem e in the mesh, similarly to point_value.

Definition at line 2110 of file system.C.

References libMesh::TypeVector< T >::add_scaled(), libMesh::FEGenericBase< Real >::build(), libMesh::Elem::contains_point(), libMesh::System::current_solution(), libMesh::Elem::dim(), libMesh::System::get_dof_map(), libMesh::FEInterface::inverse_map(), libMesh::processor_id(), and libMesh::DofObject::processor_id().

02111 {
02112   libmesh_assert_equal_to (e.processor_id(), libMesh::processor_id());
02113 
02114   // Ensuring that the given point is really in the element is an
02115   // expensive assert, but as long as debugging is turned on we might
02116   // as well try to catch a particularly nasty potential error
02117   libmesh_assert (e.contains_point(p));
02118 
02119   // Get the dof map to get the proper indices for our computation
02120   const DofMap& dof_map = this->get_dof_map();
02121 
02122   // Need dof_indices for phi[i][j]
02123   std::vector<dof_id_type> dof_indices;
02124 
02125   // Fill in the dof_indices for our element
02126   dof_map.dof_indices (&e, dof_indices, var);
02127 
02128   // Get the no of dofs assciated with this point
02129   const unsigned int num_dofs = libmesh_cast_int<unsigned int>
02130     (dof_indices.size());
02131 
02132   FEType fe_type = dof_map.variable_type(0);
02133 
02134   // Build a FE again so we can calculate u(p)
02135   AutoPtr<FEBase> fe (FEBase::build(e.dim(), fe_type));
02136 
02137   // Map the physical co-ordinates to the master co-ordinates using the inverse_map from fe_interface.h
02138   // Build a vector of point co-ordinates to send to reinit
02139   std::vector<Point> coor(1, FEInterface::inverse_map(e.dim(), fe_type, &e, p));
02140 
02141   // Get the values of the shape function derivatives
02142   const std::vector<std::vector<RealGradient> >&  dphi = fe->get_dphi();
02143 
02144   // Reinitialize the element and compute the shape function values at coor
02145   fe->reinit (&e, &coor);
02146 
02147   // Get ready to accumulate a gradient
02148   Gradient grad_u;
02149 
02150   for (unsigned int l=0; l<num_dofs; l++)
02151     {
02152       grad_u.add_scaled (dphi[l][0], this->current_solution (dof_indices[l]));
02153     }
02154 
02155   return grad_u;
02156 }

Gradient libMesh::System::point_gradient	(	unsigned int	var,
		const Point &	p,
		const bool	insist_on_success = `true`
	)			const `[inherited]`

Returns the gradient of the solution variable var at the physical point p in the mesh, similarly to point_value.

Definition at line 2062 of file system.C.

References libMesh::Parallel::Communicator::broadcast(), libMesh::CommWorld, libMesh::PointLocatorBase::enable_out_of_mesh_mode(), libMesh::System::get_mesh(), mesh, libMesh::Parallel::Communicator::min(), libMesh::n_processors(), libMesh::processor_id(), libMesh::DofObject::processor_id(), libMesh::MeshBase::sub_point_locator(), and libMesh::Parallel::Communicator::verify().

02063 {
02064   // This function must be called on every processor; there's no
02065   // telling where in the partition p falls.
02066   parallel_only();
02067 
02068   // And every processor had better agree about which point we're
02069   // looking for
02070 #ifndef NDEBUG
02071   CommWorld.verify(p);
02072 #endif // NDEBUG
02073 
02074   // Get a reference to the mesh object associated with the system object that calls this function
02075   const MeshBase &mesh = this->get_mesh();
02076 
02077   // Use an existing PointLocator or create a new one
02078   AutoPtr<PointLocatorBase> locator_ptr = mesh.sub_point_locator();
02079   PointLocatorBase& locator = *locator_ptr;
02080 
02081   if (!insist_on_success)
02082     locator.enable_out_of_mesh_mode();
02083 
02084   // Get a pointer to the element that contains P
02085   const Elem *e = locator(p);
02086 
02087   Gradient grad_u;
02088 
02089   if (e && e->processor_id() == libMesh::processor_id())
02090     grad_u = point_gradient(var, p, *e);
02091 
02092   // If I have an element containing p, then let's let everyone know
02093   processor_id_type lowest_owner =
02094     (e && (e->processor_id() == libMesh::processor_id())) ?
02095     libMesh::processor_id() : libMesh::n_processors();
02096   CommWorld.min(lowest_owner);
02097 
02098   // Everybody should get their value from a processor that was able
02099   // to compute it.
02100   // If nobody admits owning the point, we may have a problem.
02101   if (lowest_owner != libMesh::n_processors())
02102     CommWorld.broadcast(grad_u, lowest_owner);
02103   else
02104     libmesh_assert(!insist_on_success);
02105 
02106   return grad_u;
02107 }

Tensor libMesh::System::point_hessian	(	unsigned int	var,
		const Point &	p,
		const Elem &	e
	)			const `[inherited]`

Returns the second derivative tensor of the solution variable var at the physical point p in local Elem e in the mesh, similarly to point_value.

Definition at line 2208 of file system.C.

References libMesh::TypeTensor< T >::add_scaled(), libMesh::FEGenericBase< Real >::build(), libMesh::Elem::contains_point(), libMesh::System::current_solution(), libMesh::Elem::dim(), libMesh::System::get_dof_map(), libMesh::FEInterface::inverse_map(), libMesh::processor_id(), and libMesh::DofObject::processor_id().

02209 {
02210   libmesh_assert_equal_to (e.processor_id(), libMesh::processor_id());
02211 
02212   // Ensuring that the given point is really in the element is an
02213   // expensive assert, but as long as debugging is turned on we might
02214   // as well try to catch a particularly nasty potential error
02215   libmesh_assert (e.contains_point(p));
02216 
02217   // Get the dof map to get the proper indices for our computation
02218   const DofMap& dof_map = this->get_dof_map();
02219 
02220   // Need dof_indices for phi[i][j]
02221   std::vector<dof_id_type> dof_indices;
02222 
02223   // Fill in the dof_indices for our element
02224   dof_map.dof_indices (&e, dof_indices, var);
02225 
02226   // Get the no of dofs assciated with this point
02227   const unsigned int num_dofs = libmesh_cast_int<unsigned int>
02228     (dof_indices.size());
02229 
02230   FEType fe_type = dof_map.variable_type(0);
02231 
02232   // Build a FE again so we can calculate u(p)
02233   AutoPtr<FEBase> fe (FEBase::build(e.dim(), fe_type));
02234 
02235   // Map the physical co-ordinates to the master co-ordinates using the inverse_map from fe_interface.h
02236   // Build a vector of point co-ordinates to send to reinit
02237   std::vector<Point> coor(1, FEInterface::inverse_map(e.dim(), fe_type, &e, p));
02238 
02239   // Get the values of the shape function derivatives
02240   const std::vector<std::vector<RealTensor> >&  d2phi = fe->get_d2phi();
02241 
02242   // Reinitialize the element and compute the shape function values at coor
02243   fe->reinit (&e, &coor);
02244 
02245   // Get ready to accumulate a hessian
02246   Tensor hess_u;
02247 
02248   for (unsigned int l=0; l<num_dofs; l++)
02249     {
02250       hess_u.add_scaled (d2phi[l][0], this->current_solution (dof_indices[l]));
02251     }
02252 
02253   return hess_u;
02254 }

Tensor libMesh::System::point_hessian	(	unsigned int	var,
		const Point &	p,
		const bool	insist_on_success = `true`
	)			const `[inherited]`

Returns the second derivative tensor of the solution variable var at the physical point p in the mesh, similarly to point_value.

Definition at line 2161 of file system.C.

References libMesh::Parallel::Communicator::broadcast(), libMesh::CommWorld, libMesh::PointLocatorBase::enable_out_of_mesh_mode(), libMesh::System::get_mesh(), mesh, libMesh::Parallel::Communicator::min(), libMesh::n_processors(), libMesh::processor_id(), libMesh::DofObject::processor_id(), libMesh::MeshBase::sub_point_locator(), and libMesh::Parallel::Communicator::verify().

02162 {
02163   // This function must be called on every processor; there's no
02164   // telling where in the partition p falls.
02165   parallel_only();
02166 
02167   // And every processor had better agree about which point we're
02168   // looking for
02169 #ifndef NDEBUG
02170   CommWorld.verify(p);
02171 #endif // NDEBUG
02172 
02173   // Get a reference to the mesh object associated with the system object that calls this function
02174   const MeshBase &mesh = this->get_mesh();
02175 
02176   // Use an existing PointLocator or create a new one
02177   AutoPtr<PointLocatorBase> locator_ptr = mesh.sub_point_locator();
02178   PointLocatorBase& locator = *locator_ptr;
02179 
02180   if (!insist_on_success)
02181     locator.enable_out_of_mesh_mode();
02182 
02183   // Get a pointer to the element that contains P
02184   const Elem *e = locator(p);
02185 
02186   Tensor hess_u;
02187 
02188   if (e && e->processor_id() == libMesh::processor_id())
02189     hess_u = point_hessian(var, p, *e);
02190 
02191   // If I have an element containing p, then let's let everyone know
02192   processor_id_type lowest_owner =
02193     (e && (e->processor_id() == libMesh::processor_id())) ?
02194     libMesh::processor_id() : libMesh::n_processors();
02195   CommWorld.min(lowest_owner);
02196 
02197   // Everybody should get their value from a processor that was able
02198   // to compute it.
02199   // If nobody admits owning the point, we may have a problem.
02200   if (lowest_owner != libMesh::n_processors())
02201     CommWorld.broadcast(hess_u, lowest_owner);
02202   else
02203     libmesh_assert(!insist_on_success);
02204 
02205   return hess_u;
02206 }

Number libMesh::System::point_value	(	unsigned int	var,
		const Point &	p,
		const Elem &	e
	)			const `[inherited]`

Returns the value of the solution variable var at the physical point p contained in local Elem e

This version of point_value can be run in serial, but assumes e is in the local mesh partition.

Definition at line 2012 of file system.C.

References libMesh::FEGenericBase< Real >::build(), libMesh::Elem::contains_point(), libMesh::System::current_solution(), libMesh::Elem::dim(), libMesh::System::get_dof_map(), libMesh::FEInterface::inverse_map(), libMesh::processor_id(), and libMesh::DofObject::processor_id().

02013 {
02014   libmesh_assert_equal_to (e.processor_id(), libMesh::processor_id());
02015 
02016   // Ensuring that the given point is really in the element is an
02017   // expensive assert, but as long as debugging is turned on we might
02018   // as well try to catch a particularly nasty potential error
02019   libmesh_assert (e.contains_point(p));
02020 
02021   // Get the dof map to get the proper indices for our computation
02022   const DofMap& dof_map = this->get_dof_map();
02023 
02024   // Need dof_indices for phi[i][j]
02025   std::vector<dof_id_type> dof_indices;
02026 
02027   // Fill in the dof_indices for our element
02028   dof_map.dof_indices (&e, dof_indices, var);
02029 
02030   // Get the no of dofs assciated with this point
02031   const unsigned int num_dofs = libmesh_cast_int<unsigned int>
02032     (dof_indices.size());
02033 
02034   FEType fe_type = dof_map.variable_type(0);
02035 
02036   // Build a FE so we can calculate u(p)
02037   AutoPtr<FEBase> fe (FEBase::build(e.dim(), fe_type));
02038 
02039   // Map the physical co-ordinates to the master co-ordinates using the inverse_map from fe_interface.h
02040   // Build a vector of point co-ordinates to send to reinit
02041   std::vector<Point> coor(1, FEInterface::inverse_map(e.dim(), fe_type, &e, p));
02042 
02043   // Get the shape function values
02044   const std::vector<std::vector<Real> >& phi = fe->get_phi();
02045 
02046   // Reinitialize the element and compute the shape function values at coor
02047   fe->reinit (&e, &coor);
02048 
02049   // Get ready to accumulate a value
02050   Number u = 0;
02051 
02052   for (unsigned int l=0; l<num_dofs; l++)
02053     {
02054       u += phi[l][0]*this->current_solution (dof_indices[l]);
02055     }
02056 
02057   return u;
02058 }

Number libMesh::System::point_value	(	unsigned int	var,
		const Point &	p,
		const bool	insist_on_success = `true`
	)			const `[inherited]`

Returns the value of the solution variable var at the physical point p in the mesh, without knowing a priori which element contains p.

Note that this function uses MeshBase::sub_point_locator(); users may or may not want to call MeshBase::clear_point_locator() afterward. Also, point_locator() is expensive (N log N for initial construction, log N for evaluations). Avoid using this function in any context where you are already looping over elements.

Because the element containing p may lie on any processor, this function is parallel-only.

By default this method expects the point to reside inside the domain and will abort if no element can be found which contains . The optional parameter insist_on_success can be set to false to allow the method to return 0 when the point is not located.

Definition at line 1965 of file system.C.

References libMesh::Parallel::Communicator::broadcast(), libMesh::CommWorld, libMesh::PointLocatorBase::enable_out_of_mesh_mode(), libMesh::System::get_mesh(), mesh, libMesh::Parallel::Communicator::min(), libMesh::n_processors(), libMesh::processor_id(), libMesh::DofObject::processor_id(), libMesh::MeshBase::sub_point_locator(), and libMesh::Parallel::Communicator::verify().

01966 {
01967   // This function must be called on every processor; there's no
01968   // telling where in the partition p falls.
01969   parallel_only();
01970 
01971   // And every processor had better agree about which point we're
01972   // looking for
01973 #ifndef NDEBUG
01974   CommWorld.verify(p);
01975 #endif // NDEBUG
01976 
01977   // Get a reference to the mesh object associated with the system object that calls this function
01978   const MeshBase &mesh = this->get_mesh();
01979 
01980   // Use an existing PointLocator or create a new one
01981   AutoPtr<PointLocatorBase> locator_ptr = mesh.sub_point_locator();
01982   PointLocatorBase& locator = *locator_ptr;
01983 
01984   if (!insist_on_success)
01985     locator.enable_out_of_mesh_mode();
01986 
01987   // Get a pointer to the element that contains P
01988   const Elem *e = locator(p);
01989 
01990   Number u = 0;
01991 
01992   if (e && e->processor_id() == libMesh::processor_id())
01993     u = point_value(var, p, *e);
01994 
01995   // If I have an element containing p, then let's let everyone know
01996   processor_id_type lowest_owner =
01997     (e && (e->processor_id() == libMesh::processor_id())) ?
01998     libMesh::processor_id() : libMesh::n_processors();
01999   CommWorld.min(lowest_owner);
02000 
02001   // Everybody should get their value from a processor that was able
02002   // to compute it.
02003   // If nobody admits owning the point, we have a problem.
02004   if (lowest_owner != libMesh::n_processors())
02005     CommWorld.broadcast(u, lowest_owner);
02006   else
02007     libmesh_assert(!insist_on_success);
02008 
02009   return u;
02010 }

void libMesh::ReferenceCounter::print_info ( std::ostream & out = libMesh::out ) [static, inherited]

Prints the reference information, by default to libMesh::out.

Definition at line 88 of file reference_counter.C.

References libMesh::ReferenceCounter::_enable_print_counter, and libMesh::ReferenceCounter::get_info().

00089 {
00090   if( _enable_print_counter ) out_stream << ReferenceCounter::get_info();
00091 }

void libMesh::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,
		const Parameters &	parameters
	)			const `[inherited]`

Projects arbitrary functions onto the current solution. The function value fptr and its gradient gptr are represented by function pointers. A gradient gptr is only required/used for projecting onto finite element spaces with continuous derivatives.

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

Definition at line 450 of file system_projection.C.

References libMesh::System::project_solution().

00459 {
00460   WrappedFunction<Number> f(*this, fptr, &parameters);
00461   WrappedFunction<Gradient> g(*this, gptr, &parameters);
00462   this->project_solution(&f, &g);
00463 }

void libMesh::System::project_solution	(	FEMFunctionBase< Number > *	f,
		FEMFunctionBase< Gradient > *	g = `NULL`
	)			const `[inherited]`

Projects arbitrary functions onto the current solution. The function value f and its gradient g are user-provided cloneable functors. A gradient g is only required/used for projecting onto finite element spaces with continuous derivatives. If non-default Parameters are to be used, they can be provided in the parameters argument.

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

Definition at line 483 of file system_projection.C.

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

00485 {
00486   this->project_vector(*solution, f, g);
00487 
00488   solution->localize(*current_local_solution, _dof_map->get_send_list());
00489 }

void libMesh::System::project_solution	(	FunctionBase< Number > *	f,
		FunctionBase< Gradient > *	g = `NULL`
	)			const `[inherited]`

Projects arbitrary functions onto the current solution. The function value f and its gradient g are user-provided cloneable functors. A gradient g is only required/used for projecting onto finite element spaces with continuous derivatives. If non-default Parameters are to be used, they can be provided in the parameters argument.

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

Definition at line 470 of file system_projection.C.

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

Referenced by libMesh::System::project_solution().

00472 {
00473   this->project_vector(*solution, f, g);
00474 
00475   solution->localize(*current_local_solution, _dof_map->get_send_list());
00476 }

bool& libMesh::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 761 of file system.h.

References libMesh::System::_solution_projection.

Referenced by libMesh::AdjointRefinementEstimator::estimate_error(), and libMesh::MemorySolutionHistory::store().

00762     { return _solution_projection; }

void libMesh::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 272 of file system_projection.C.

00274 {
00275   START_LOG ("project_vector()", "System");
00276 
00283   new_v.clear();
00284 
00285 #ifdef LIBMESH_ENABLE_AMR
00286 
00287   // Resize the new vector and get a serial version.
00288   NumericVector<Number> *new_vector_ptr = NULL;
00289   AutoPtr<NumericVector<Number> > new_vector_built;
00290   NumericVector<Number> *local_old_vector;
00291   AutoPtr<NumericVector<Number> > local_old_vector_built;
00292   const NumericVector<Number> *old_vector_ptr = NULL;
00293 
00294   ConstElemRange active_local_elem_range
00295     (this->get_mesh().active_local_elements_begin(),
00296      this->get_mesh().active_local_elements_end());
00297 
00298   // If the old vector was uniprocessor, make the new
00299   // vector uniprocessor
00300   if (old_v.type() == SERIAL)
00301     {
00302       new_v.init (this->n_dofs(), false, SERIAL);
00303       new_vector_ptr = &new_v;
00304       old_vector_ptr = &old_v;
00305     }
00306 
00307   // Otherwise it is a parallel, distributed vector, which
00308   // we need to localize.
00309   else if (old_v.type() == PARALLEL)
00310     {
00311       // Build a send list for efficient localization
00312       BuildProjectionList projection_list(*this);
00313       Threads::parallel_reduce (active_local_elem_range,
00314                                 projection_list);
00315 
00316       // Create a sorted, unique send_list
00317       projection_list.unique();
00318 
00319       new_v.init (this->n_dofs(), this->n_local_dofs(), false, PARALLEL);
00320       new_vector_built = NumericVector<Number>::build();
00321       local_old_vector_built = NumericVector<Number>::build();
00322       new_vector_ptr = new_vector_built.get();
00323       local_old_vector = local_old_vector_built.get();
00324       new_vector_ptr->init(this->n_dofs(), false, SERIAL);
00325       local_old_vector->init(old_v.size(), false, SERIAL);
00326       old_v.localize(*local_old_vector, projection_list.send_list);
00327       local_old_vector->close();
00328       old_vector_ptr = local_old_vector;
00329     }
00330   else if (old_v.type() == GHOSTED)
00331     {
00332       // Build a send list for efficient localization
00333       BuildProjectionList projection_list(*this);
00334       Threads::parallel_reduce (active_local_elem_range,
00335                                 projection_list);
00336 
00337       // Create a sorted, unique send_list
00338       projection_list.unique();
00339 
00340       new_v.init (this->n_dofs(), this->n_local_dofs(),
00341                   this->get_dof_map().get_send_list(), false, GHOSTED);
00342 
00343       local_old_vector_built = NumericVector<Number>::build();
00344       new_vector_ptr = &new_v;
00345       local_old_vector = local_old_vector_built.get();
00346       local_old_vector->init(old_v.size(), old_v.local_size(),
00347                              projection_list.send_list, false, GHOSTED);
00348       old_v.localize(*local_old_vector, projection_list.send_list);
00349       local_old_vector->close();
00350       old_vector_ptr = local_old_vector;
00351     }
00352   else // unknown old_v.type()
00353     {
00354       libMesh::err << "ERROR: Unknown old_v.type() == " << old_v.type()
00355                     << std::endl;
00356       libmesh_error();
00357     }
00358 
00359   // Note that the above will have zeroed the new_vector.
00360   // Just to be sure, assert that new_vector_ptr and old_vector_ptr
00361   // were successfully set before trying to deref them.
00362   libmesh_assert(new_vector_ptr);
00363   libmesh_assert(old_vector_ptr);
00364 
00365   NumericVector<Number> &new_vector = *new_vector_ptr;
00366   const NumericVector<Number> &old_vector = *old_vector_ptr;
00367 
00368   Threads::parallel_for (active_local_elem_range,
00369                          ProjectVector(*this,
00370                                        old_vector,
00371                                        new_vector)
00372                          );
00373 
00374   // Copy the SCALAR dofs from old_vector to new_vector
00375   // Note: We assume that all SCALAR dofs are on the
00376   // processor with highest ID
00377   if(libMesh::processor_id() == (libMesh::n_processors()-1))
00378   {
00379     const DofMap& dof_map = this->get_dof_map();
00380     for (unsigned int var=0; var<this->n_vars(); var++)
00381       if(this->variable(var).type().family == SCALAR)
00382         {
00383           // We can just map SCALAR dofs directly across
00384           std::vector<dof_id_type> new_SCALAR_indices, old_SCALAR_indices;
00385           dof_map.SCALAR_dof_indices (new_SCALAR_indices, var, false);
00386           dof_map.SCALAR_dof_indices (old_SCALAR_indices, var, true);
00387           const unsigned int new_n_dofs = 
00388             libmesh_cast_int<unsigned int>(new_SCALAR_indices.size());
00389 
00390           for (unsigned int i=0; i<new_n_dofs; i++)
00391           {
00392             new_vector.set( new_SCALAR_indices[i], old_vector(old_SCALAR_indices[i]) );
00393           }
00394         }
00395   }
00396 
00397   new_vector.close();
00398 
00399   // If the old vector was serial, we probably need to send our values
00400   // to other processors
00401   //
00402   // FIXME: I'm not sure how to make a NumericVector do that without
00403   // creating a temporary parallel vector to use localize! - RHS
00404   if (old_v.type() == SERIAL)
00405     {
00406       AutoPtr<NumericVector<Number> > dist_v = NumericVector<Number>::build();
00407       dist_v->init(this->n_dofs(), this->n_local_dofs(), false, PARALLEL);
00408       dist_v->close();
00409 
00410       for (dof_id_type i=0; i!=dist_v->size(); i++)
00411         if (new_vector(i) != 0.0)
00412           dist_v->set(i, new_vector(i));
00413 
00414       dist_v->close();
00415 
00416       dist_v->localize (new_v, this->get_dof_map().get_send_list());
00417       new_v.close();
00418     }
00419   // If the old vector was parallel, we need to update it
00420   // and free the localized copies
00421   else if (old_v.type() == PARALLEL)
00422     {
00423       // We may have to set dof values that this processor doesn't
00424       // own in certain special cases, like LAGRANGE FIRST or
00425       // HERMITE THIRD elements on second-order meshes
00426       for (dof_id_type i=0; i!=new_v.size(); i++)
00427         if (new_vector(i) != 0.0)
00428           new_v.set(i, new_vector(i));
00429       new_v.close();
00430     }
00431 
00432   this->get_dof_map().enforce_constraints_exactly(*this, &new_v);
00433 
00434 #else
00435 
00436   // AMR is disabled: simply copy the vector
00437   new_v = old_v;
00438 
00439 #endif // #ifdef LIBMESH_ENABLE_AMR
00440 
00441   STOP_LOG("project_vector()", "System");
00442 }

void libMesh::System::project_vector ( NumericVector< Number > & vector ) const [protected, inherited]

Projects the vector defined on the old mesh onto the new mesh.

Definition at line 255 of file system_projection.C.

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

00256 {
00257   // Create a copy of the vector, which currently
00258   // contains the old data.
00259   AutoPtr<NumericVector<Number> >
00260     old_vector (vector.clone());
00261 
00262   // Project the old vector to the new vector
00263   this->project_vector (*old_vector, vector);
00264 }

void libMesh::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,
		const Parameters &	parameters,
		NumericVector< Number > &	new_vector
	)			const `[inherited]`

Projects arbitrary functions onto a vector of degree of freedom values for the current system. The function value fptr and its gradient gptr are represented by function pointers. A gradient gptr is only required/used for projecting onto finite element spaces with continuous derivatives.

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

Definition at line 496 of file system_projection.C.

References libMesh::System::project_vector().

00506 {
00507   WrappedFunction<Number> f(*this, fptr, &parameters);
00508   WrappedFunction<Gradient> g(*this, gptr, &parameters);
00509   this->project_vector(new_vector, &f, &g);
00510 }

void libMesh::System::project_vector	(	NumericVector< Number > &	new_vector,
		FEMFunctionBase< Number > *	f,
		FEMFunctionBase< Gradient > *	g = `NULL`
	)			const `[inherited]`

Projects arbitrary functions onto a vector of degree of freedom values for the current system. The function value f and its gradient g are user-provided cloneable functors. A gradient g is only required/used for projecting onto finite element spaces with continuous derivatives. If non-default Parameters are to be used, they can be provided in the parameters argument.

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

Definition at line 579 of file system_projection.C.

References libMesh::MeshBase::active_local_elements_begin(), libMesh::NumericVector< T >::close(), libMesh::FEMFunctionBase< Output >::component(), libMesh::DofMap::enforce_constraints_exactly(), libMesh::FEType::family, libMesh::System::get_dof_map(), libMesh::System::get_mesh(), libMesh::n_processors(), libMesh::System::n_vars(), libMesh::Threads::parallel_for(), libMesh::FEMContext::pre_fe_reinit(), libMesh::processor_id(), libMeshEnums::SCALAR, libMesh::DofMap::SCALAR_dof_indices(), libMesh::NumericVector< T >::set(), libMesh::System::time, libMesh::Variable::type(), libMesh::System::variable(), and libMesh::System::variable_scalar_number().

00582 {
00583   START_LOG ("project_fem_vector()", "System");
00584 
00585   Threads::parallel_for
00586     (ConstElemRange (this->get_mesh().active_local_elements_begin(),
00587                      this->get_mesh().active_local_elements_end() ),
00588      ProjectFEMSolution(*this, f, g, new_vector)
00589     );
00590 
00591   // Also, load values into the SCALAR dofs
00592   // Note: We assume that all SCALAR dofs are on the
00593   // processor with highest ID
00594   if(libMesh::processor_id() == (libMesh::n_processors()-1))
00595   {
00596     // FIXME: Do we want to first check for SCALAR vars before building this? [PB]
00597     FEMContext context( *this );
00598 
00599     const DofMap& dof_map = this->get_dof_map();
00600     for (unsigned int var=0; var<this->n_vars(); var++)
00601       if(this->variable(var).type().family == SCALAR)
00602         {
00603           // FIXME: We reinit with an arbitrary element in case the user
00604           //        doesn't override FEMFunctionBase::component. Is there
00605           //        any use case we're missing? [PB]
00606           Elem *el = const_cast<Elem *>(*(this->get_mesh().active_local_elements_begin()));
00607           context.pre_fe_reinit( *this, el );
00608 
00609           std::vector<dof_id_type> SCALAR_indices;
00610           dof_map.SCALAR_dof_indices (SCALAR_indices, var);
00611           const unsigned int n_SCALAR_dofs =
00612             libmesh_cast_int<unsigned int>(SCALAR_indices.size());
00613 
00614           for (unsigned int i=0; i<n_SCALAR_dofs; i++)
00615           {
00616             const dof_id_type global_index = SCALAR_indices[i];
00617             const unsigned int component_index =
00618               this->variable_scalar_number(var,i);
00619             
00620             new_vector.set(global_index, f->component(context, component_index, Point(), this->time));
00621           }
00622         }
00623   }
00624 
00625   new_vector.close();
00626 
00627 #ifdef LIBMESH_ENABLE_CONSTRAINTS
00628   this->get_dof_map().enforce_constraints_exactly(*this, &new_vector);
00629 #endif
00630 
00631   STOP_LOG("project_fem_vector()", "System");
00632 }

void libMesh::System::project_vector	(	NumericVector< Number > &	new_vector,
		FunctionBase< Number > *	f,
		FunctionBase< Gradient > *	g = `NULL`
	)			const `[inherited]`

Projects arbitrary functions onto a vector of degree of freedom values for the current system. The function value f and its gradient g are user-provided cloneable functors. A gradient g is only required/used for projecting onto finite element spaces with continuous derivatives. If non-default Parameters are to be used, they can be provided in the parameters argument.

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

Definition at line 516 of file system_projection.C.

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

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

00519 {
00520   START_LOG ("project_vector()", "System");
00521 
00522   Threads::parallel_for
00523     (ConstElemRange (this->get_mesh().active_local_elements_begin(),
00524                      this->get_mesh().active_local_elements_end() ),
00525      ProjectSolution(*this, f, g,
00526                      this->get_equation_systems().parameters,
00527                      new_vector)
00528     );
00529 
00530   // Also, load values into the SCALAR dofs
00531   // Note: We assume that all SCALAR dofs are on the
00532   // processor with highest ID
00533   if(libMesh::processor_id() == (libMesh::n_processors()-1))
00534   {
00535     // We get different scalars as different
00536     // components from a new-style f functor.
00537     DenseVector<Number> fout(this->n_components());
00538     bool filled_fout = false;
00539 
00540     const DofMap& dof_map = this->get_dof_map();
00541     for (unsigned int var=0; var<this->n_vars(); var++)
00542       if(this->variable(var).type().family == SCALAR)
00543         {
00544           if (!filled_fout)
00545             {
00546               (*f) (Point(), this->time, fout);
00547               filled_fout = true;
00548             }
00549 
00550           std::vector<dof_id_type> SCALAR_indices;
00551           dof_map.SCALAR_dof_indices (SCALAR_indices, var);
00552           const unsigned int n_SCALAR_dofs = 
00553             libmesh_cast_int<unsigned int>(SCALAR_indices.size());
00554 
00555           for (unsigned int i=0; i<n_SCALAR_dofs; i++)
00556           {
00557             const dof_id_type global_index = SCALAR_indices[i];
00558             const unsigned int component_index =
00559               this->variable_scalar_number(var,i);
00560             new_vector.set(global_index, fout(component_index));
00561           }
00562         }
00563   }
00564 
00565   new_vector.close();
00566 
00567 #ifdef LIBMESH_ENABLE_CONSTRAINTS
00568   this->get_dof_map().enforce_constraints_exactly(*this, &new_vector);
00569 #endif
00570 
00571   STOP_LOG("project_vector()", "System");
00572 }

void libMesh::System::prolong_vectors ( ) [virtual, inherited]

Prolong vectors after the mesh has refined

Definition at line 368 of file system.C.

References libMesh::System::restrict_vectors().

Referenced by libMesh::EquationSystems::reinit().

00369 {
00370 #ifdef LIBMESH_ENABLE_AMR
00371   // Currently project_vector handles both restriction and prolongation
00372   this->restrict_vectors();
00373 #endif
00374 }

void libMesh::ImplicitSystem::qoi_parameter_hessian	(	const QoISet &	qoi_indices,
		const ParameterVector &	parameters,
		SensitivityData &	hessian
	)			`[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) This Hessian is the output of this method, where for each q_i, H_jk is stored in hessian.second_derivative(i,j,k).

Reimplemented from libMesh::System.

Definition at line 1080 of file implicit_system.C.

References libMesh::SensitivityData::allocate_hessian_data(), libMesh::QoISet::has_index(), libMesh::Real, libMesh::SensitivityData::second_derivative(), libMesh::ParameterVector::size(), and libMesh::TOLERANCE.

01083 {
01084   // We currently get partial derivatives via finite differencing
01085   const Real delta_p = TOLERANCE;
01086 
01087   // We'll use one temporary vector for matrix-vector-vector products
01088   AutoPtr<NumericVector<Number> > tempvec = this->solution->zero_clone();
01089 
01090   // And another temporary vector to hold a copy of the true solution
01091   // so we can safely perturb this->solution.
01092   AutoPtr<NumericVector<Number> > oldsolution = this->solution->clone();
01093 
01094   const unsigned int Np = libmesh_cast_int<unsigned int>
01095     (parameters.size());
01096   const unsigned int Nq = libmesh_cast_int<unsigned int>
01097     (qoi.size());
01098 
01099   // For each quantity of interest q, the parameter sensitivity
01100   // Hessian is defined as q''_{kl} = {d^2 q}/{d p_k d p_l}.
01101   //
01102   // We calculate it from values and partial derivatives of the
01103   // quantity of interest function Q, solution u, adjoint solution z,
01104   // and residual R, as:
01105   //
01106   // q''_{kl} =
01107   // Q''_{kl} + Q''_{uk}(u)*u'_l + Q''_{ul}(u) * u'_k +
01108   // Q''_{uu}(u)*u'_k*u'_l -
01109   // R''_{kl}(u,z) -
01110   // R''_{uk}(u,z)*u'_l - R''_{ul}(u,z)*u'_k -
01111   // R''_{uu}(u,z)*u'_k*u'_l
01112   //
01113   // See the adjoints model document for more details.
01114 
01115   // We first do an adjoint solve to get z for each quantity of
01116   // interest
01117   // if we havent already or dont have an initial condition for the adjoint
01118   if (!this->is_adjoint_already_solved())
01119     {
01120       this->adjoint_solve(qoi_indices);
01121     }
01122 
01123   // And a sensitivity solve to get u_k for each parameter
01124   this->sensitivity_solve(parameters);
01125 
01126   // Get ready to fill in second derivatives:
01127   sensitivities.allocate_hessian_data(qoi_indices, *this, parameters);
01128 
01129   for (unsigned int k=0; k != Np; ++k)
01130     {
01131       Number old_parameterk = *parameters[k];
01132 
01133       // The Hessian is symmetric, so we just calculate the lower
01134       // triangle and the diagonal, and we get the upper triangle from
01135       // the transpose of the lower
01136 
01137       for (unsigned int l=0; l != k+1; ++l)
01138         {
01139           // The second partial derivatives with respect to parameters
01140           // are all calculated via a central finite difference
01141           // stencil:
01142           // F''_{kl} ~= (F(p+dp*e_k+dp*e_l) - F(p+dp*e_k-dp*e_l) -
01143           //              F(p-dp*e_k+dp*e_l) + F(p-dp*e_k-dp*e_l))/(4*dp^2)
01144           // We will add Q''_{kl}(u) and subtract R''_{kl}(u,z) at the
01145           // same time.
01146           //
01147           // We have to be careful with the perturbations to handle
01148           // the k=l case
01149 
01150           Number old_parameterl = *parameters[l];
01151 
01152           *parameters[k] += delta_p;
01153           *parameters[l] += delta_p;
01154           this->assemble_qoi(qoi_indices);
01155           this->assembly(true, false);
01156           this->rhs->close();
01157           std::vector<Number> partial2q_term = this->qoi;
01158           std::vector<Number> partial2R_term(this->qoi.size());
01159           for (unsigned int i=0; i != Nq; ++i)
01160             if (qoi_indices.has_index(i))
01161               partial2R_term[i] = this->rhs->dot(this->get_adjoint_solution(i));
01162 
01163           *parameters[l] -= 2.*delta_p;
01164           this->assemble_qoi(qoi_indices);
01165           this->assembly(true, false);
01166           this->rhs->close();
01167           for (unsigned int i=0; i != Nq; ++i)
01168             if (qoi_indices.has_index(i))
01169               {
01170                 partial2q_term[i] -= this->qoi[i];
01171                 partial2R_term[i] -= this->rhs->dot(this->get_adjoint_solution(i));
01172               }
01173 
01174           *parameters[k] -= 2.*delta_p;
01175           this->assemble_qoi(qoi_indices);
01176           this->assembly(true, false);
01177           this->rhs->close();
01178           for (unsigned int i=0; i != Nq; ++i)
01179             if (qoi_indices.has_index(i))
01180               {
01181                 partial2q_term[i] += this->qoi[i];
01182                 partial2R_term[i] += this->rhs->dot(this->get_adjoint_solution(i));
01183               }
01184 
01185           *parameters[l] += 2.*delta_p;
01186           this->assemble_qoi(qoi_indices);
01187           this->assembly(true, false);
01188           this->rhs->close();
01189           for (unsigned int i=0; i != Nq; ++i)
01190             if (qoi_indices.has_index(i))
01191               {
01192                 partial2q_term[i] -= this->qoi[i];
01193                 partial2R_term[i] -= this->rhs->dot(this->get_adjoint_solution(i));
01194                 partial2q_term[i] /= (4. * delta_p * delta_p);
01195                 partial2R_term[i] /= (4. * delta_p * delta_p);
01196               }
01197 
01198           for (unsigned int i=0; i != Nq; ++i)
01199             if (qoi_indices.has_index(i))
01200               {
01201                 Number current_terms = partial2q_term[i] - partial2R_term[i];
01202                 sensitivities.second_derivative(i,k,l) += current_terms;
01203                 if (k != l)
01204                   sensitivities.second_derivative(i,l,k) += current_terms;
01205               }
01206 
01207           // Don't leave the parameters perturbed
01208           *parameters[l] = old_parameterl;
01209           *parameters[k] = old_parameterk;
01210         }
01211 
01212       // We get (partial q / partial u) and
01213       // (partial R / partial u) from the user, but centrally
01214       // difference to get q_uk and R_uk terms:
01215       // (partial^2 q / partial u partial k)
01216       // q_uk*u'_l = (q_u(p+dp*e_k)*u'_l - q_u(p-dp*e_k)*u'_l)/(2*dp)
01217       // R_uk*z*u'_l = (R_u(p+dp*e_k)*z*u'_l - R_u(p-dp*e_k)*z*u'_l)/(2*dp)
01218       //
01219       // To avoid creating Nq temporary vectors, we add these
01220       // subterms to the sensitivities output one by one.
01221       //
01222       // FIXME: this is probably a bad order of operations for
01223       // controlling floating point error.
01224 
01225       *parameters[k] = old_parameterk + delta_p;
01226       this->assembly(false, true);
01227       this->matrix->close();
01228       this->assemble_qoi_derivative(qoi_indices);
01229 
01230       for (unsigned int l=0; l != Np; ++l)
01231         {
01232           this->matrix->vector_mult(*tempvec, this->get_sensitivity_solution(l));
01233           for (unsigned int i=0; i != Nq; ++i)
01234             if (qoi_indices.has_index(i))
01235               {
01236                 this->get_adjoint_rhs(i).close();
01237                 Number current_terms =
01238                   (this->get_adjoint_rhs(i).dot(this->get_sensitivity_solution(l)) -
01239                    tempvec->dot(this->get_adjoint_solution(i))) / (2.*delta_p);
01240                 sensitivities.second_derivative(i,k,l) += current_terms;
01241 
01242                 // We use the _uk terms twice; symmetry lets us reuse
01243                 // these calculations for the _ul terms.
01244 
01245                 sensitivities.second_derivative(i,l,k) += current_terms;
01246               }
01247         }
01248 
01249       *parameters[k] = old_parameterk - delta_p;
01250       this->assembly(false, true);
01251       this->matrix->close();
01252       this->assemble_qoi_derivative(qoi_indices);
01253 
01254       for (unsigned int l=0; l != Np; ++l)
01255         {
01256           this->matrix->vector_mult(*tempvec, this->get_sensitivity_solution(l));
01257           for (unsigned int i=0; i != Nq; ++i)
01258             if (qoi_indices.has_index(i))
01259               {
01260                 this->get_adjoint_rhs(i).close();
01261                 Number current_terms =
01262                   (-this->get_adjoint_rhs(i).dot(this->get_sensitivity_solution(l)) +
01263                    tempvec->dot(this->get_adjoint_solution(i))) / (2.*delta_p);
01264                 sensitivities.second_derivative(i,k,l) += current_terms;
01265 
01266                 // We use the _uk terms twice; symmetry lets us reuse
01267                 // these calculations for the _ul terms.
01268 
01269                 sensitivities.second_derivative(i,l,k) += current_terms;
01270               }
01271         }
01272 
01273       // Don't leave the parameter perturbed
01274       *parameters[k] = old_parameterk;
01275 
01276       // Our last remaining terms are -R_uu(u,z)*u_k*u_l and
01277       // Q_uu(u)*u_k*u_l
01278       //
01279       // We take directional central finite differences of R_u and Q_u
01280       // to approximate these terms, e.g.:
01281       //
01282       // Q_uu(u)*u_k ~= (Q_u(u+dp*u_k) - Q_u(u-dp*u_k))/(2*dp)
01283 
01284       *this->solution = this->get_sensitivity_solution(k);
01285       *this->solution *= delta_p;
01286       *this->solution += *oldsolution;
01287       this->assembly(false, true);
01288       this->matrix->close();
01289       this->assemble_qoi_derivative(qoi_indices);
01290 
01291       // The Hessian is symmetric, so we just calculate the lower
01292       // triangle and the diagonal, and we get the upper triangle from
01293       // the transpose of the lower
01294       //
01295       // Note that, because we took the directional finite difference
01296       // with respect to k and not l, we've added an O(delta_p^2)
01297       // error to any permutational symmetry in the Hessian...
01298       for (unsigned int l=0; l != k+1; ++l)
01299         {
01300           this->matrix->vector_mult(*tempvec, this->get_sensitivity_solution(l));
01301           for (unsigned int i=0; i != Nq; ++i)
01302             if (qoi_indices.has_index(i))
01303               {
01304                 this->get_adjoint_rhs(i).close();
01305                 Number current_terms =
01306                   (this->get_adjoint_rhs(i).dot(this->get_sensitivity_solution(l)) -
01307                    tempvec->dot(this->get_adjoint_solution(i))) / (2.*delta_p);
01308                 sensitivities.second_derivative(i,k,l) += current_terms;
01309                 if (k != l)
01310                   sensitivities.second_derivative(i,l,k) += current_terms;
01311               }
01312         }
01313 
01314       *this->solution = this->get_sensitivity_solution(k);
01315       *this->solution *= -delta_p;
01316       *this->solution += *oldsolution;
01317       this->assembly(false, true);
01318       this->matrix->close();
01319       this->assemble_qoi_derivative(qoi_indices);
01320 
01321       for (unsigned int l=0; l != k+1; ++l)
01322         {
01323           this->matrix->vector_mult(*tempvec, this->get_sensitivity_solution(l));
01324           for (unsigned int i=0; i != Nq; ++i)
01325             if (qoi_indices.has_index(i))
01326               {
01327                 this->get_adjoint_rhs(i).close();
01328                 Number current_terms =
01329                   (-this->get_adjoint_rhs(i).dot(this->get_sensitivity_solution(l)) +
01330                    tempvec->dot(this->get_adjoint_solution(i))) / (2.*delta_p);
01331                 sensitivities.second_derivative(i,k,l) += current_terms;
01332                 if (k != l)
01333                   sensitivities.second_derivative(i,l,k) += current_terms;
01334               }
01335         }
01336 
01337       // Don't leave the solution perturbed
01338       *this->solution = *oldsolution;
01339     }
01340 
01341   // All parameters have been reset.
01342   // Don't leave the qoi or system changed - principle of least
01343   // surprise.
01344   this->assembly(true, true);
01345   this->rhs->close();
01346   this->matrix->close();
01347   this->assemble_qoi(qoi_indices);
01348 }

void libMesh::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 libMesh::System.

Definition at line 880 of file implicit_system.C.

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

00884 {
00885   // We currently get partial derivatives via finite differencing
00886   const Real delta_p = TOLERANCE;
00887 
00888   // We'll use a single temporary vector for matrix-vector-vector products
00889   AutoPtr<NumericVector<Number> > tempvec = this->solution->zero_clone();
00890 
00891   const unsigned int Np = libmesh_cast_int<unsigned int>
00892     (parameters.size());
00893   const unsigned int Nq = libmesh_cast_int<unsigned int>
00894     (qoi.size());
00895 
00896   // For each quantity of interest q, the parameter sensitivity
00897   // Hessian is defined as q''_{kl} = {d^2 q}/{d p_k d p_l}.
00898   // Given a vector of parameter perturbation weights w_l, this
00899   // function evaluates the hessian-vector product sum_l(q''_{kl}*w_l)
00900   //
00901   // We calculate it from values and partial derivatives of the
00902   // quantity of interest function Q, solution u, adjoint solution z,
00903   // parameter sensitivity adjoint solutions z^l, and residual R, as:
00904   //
00905   // sum_l(q''_{kl}*w_l) =
00906   // sum_l(w_l * Q''_{kl}) + Q''_{uk}(u)*(sum_l(w_l u'_l)) -
00907   // R'_k(u, sum_l(w_l*z^l)) - R'_{uk}(u,z)*(sum_l(w_l u'_l) -
00908   // sum_l(w_l*R''_{kl}(u,z))
00909   //
00910   // See the adjoints model document for more details.
00911 
00912   // We first do an adjoint solve to get z for each quantity of
00913   // interest
00914   // if we havent already or dont have an initial condition for the adjoint
00915   if (!this->is_adjoint_already_solved())
00916     {
00917       this->adjoint_solve(qoi_indices);
00918     }
00919 
00920   // Get ready to fill in senstivities:
00921   sensitivities.allocate_data(qoi_indices, *this, parameters);
00922 
00923   // We can't solve for all the solution sensitivities u'_l or for all
00924   // of the parameter sensitivity adjoint solutions z^l without
00925   // requiring O(Nq*Np) linear solves.  So we'll solve directly for their
00926   // weighted sum - this is just O(Nq) solves.
00927 
00928   // First solve for sum_l(w_l u'_l).
00929   this->weighted_sensitivity_solve(parameters, vector);
00930 
00931   // Then solve for sum_l(w_l z^l).
00932   this->weighted_sensitivity_adjoint_solve(parameters, vector, qoi_indices);
00933 
00934   for (unsigned int k=0; k != Np; ++k)
00935     {
00936       // We approximate sum_l(w_l * Q''_{kl}) with a central
00937       // differencing pertubation:
00938       // sum_l(w_l * Q''_{kl}) ~=
00939       // (Q(p + dp*w_l*e_l + dp*e_k) - Q(p - dp*w_l*e_l + dp*e_k) -
00940       // Q(p + dp*w_l*e_l - dp*e_k) + Q(p - dp*w_l*e_l - dp*e_k))/(4*dp^2)
00941 
00942       // The sum(w_l*R''_kl) term requires the same sort of pertubation,
00943       // and so we subtract it in at the same time:
00944       // sum_l(w_l * R''_{kl}) ~=
00945       // (R(p + dp*w_l*e_l + dp*e_k) - R(p - dp*w_l*e_l + dp*e_k) -
00946       // R(p + dp*w_l*e_l - dp*e_k) + R(p - dp*w_l*e_l - dp*e_k))/(4*dp^2)
00947 
00948       ParameterVector oldparameters, parameterperturbation;
00949       parameters.deep_copy(oldparameters);
00950       vector.deep_copy(parameterperturbation);
00951       parameterperturbation *= delta_p;
00952       parameters += parameterperturbation;
00953 
00954       Number old_parameter = *parameters[k];
00955 
00956       *parameters[k] = old_parameter + delta_p;
00957       this->assemble_qoi(qoi_indices);
00958       this->assembly(true, false);
00959       this->rhs->close();
00960       std::vector<Number> partial2q_term = this->qoi;
00961       std::vector<Number> partial2R_term(this->qoi.size());
00962       for (unsigned int i=0; i != Nq; ++i)
00963         if (qoi_indices.has_index(i))
00964           partial2R_term[i] = this->rhs->dot(this->get_adjoint_solution(i));
00965 
00966       *parameters[k] = old_parameter - delta_p;
00967       this->assemble_qoi(qoi_indices);
00968       this->assembly(true, false);
00969       this->rhs->close();
00970       for (unsigned int i=0; i != Nq; ++i)
00971         if (qoi_indices.has_index(i))
00972           {
00973             partial2q_term[i] -= this->qoi[i];
00974             partial2R_term[i] -= this->rhs->dot(this->get_adjoint_solution(i));
00975           }
00976 
00977       oldparameters.value_copy(parameters);
00978       parameterperturbation *= -1.0;
00979       parameters += parameterperturbation;
00980 
00981       // Re-center old_parameter, which may be affected by vector
00982       old_parameter = *parameters[k];
00983 
00984       *parameters[k] = old_parameter + delta_p;
00985       this->assemble_qoi(qoi_indices);
00986       this->assembly(true, false);
00987       this->rhs->close();
00988       for (unsigned int i=0; i != Nq; ++i)
00989         if (qoi_indices.has_index(i))
00990           {
00991             partial2q_term[i] -= this->qoi[i];
00992             partial2R_term[i] -= this->rhs->dot(this->get_adjoint_solution(i));
00993           }
00994 
00995       *parameters[k] = old_parameter - delta_p;
00996       this->assemble_qoi(qoi_indices);
00997       this->assembly(true, false);
00998       this->rhs->close();
00999       for (unsigned int i=0; i != Nq; ++i)
01000         if (qoi_indices.has_index(i))
01001           {
01002             partial2q_term[i] += this->qoi[i];
01003             partial2R_term[i] += this->rhs->dot(this->get_adjoint_solution(i));
01004           }
01005 
01006       for (unsigned int i=0; i != Nq; ++i)
01007         if (qoi_indices.has_index(i))
01008           {
01009             partial2q_term[i] /= (4. * delta_p * delta_p);
01010             partial2R_term[i] /= (4. * delta_p * delta_p);
01011           }
01012 
01013       for (unsigned int i=0; i != Nq; ++i)
01014         if (qoi_indices.has_index(i))
01015           sensitivities[i][k] = partial2q_term[i] - partial2R_term[i];
01016 
01017       // We get (partial q / partial u), R, and
01018       // (partial R / partial u) from the user, but centrally
01019       // difference to get q_uk, R_k, and R_uk terms:
01020       // (partial R / partial k)
01021       // 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)
01022       // (partial^2 q / partial u partial k)
01023       // q_uk = (q_u(p+dp*e_k) - q_u(p-dp*e_k))/(2*dp)
01024       // (partial^2 R / partial u partial k)
01025       // R_uk*z*sum(w_l*u'_l) = (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)
01026 
01027       // To avoid creating Nq temporary vectors for q_uk or R_uk, we add
01028       // subterms to the sensitivities output one by one.
01029       //
01030       // FIXME: this is probably a bad order of operations for
01031       // controlling floating point error.
01032 
01033       *parameters[k] = old_parameter + delta_p;
01034       this->assembly(true, true);
01035       this->rhs->close();
01036       this->matrix->close();
01037       this->assemble_qoi_derivative(qoi_indices);
01038 
01039       this->matrix->vector_mult(*tempvec, this->get_weighted_sensitivity_solution());
01040 
01041       for (unsigned int i=0; i != Nq; ++i)
01042         if (qoi_indices.has_index(i))
01043           {
01044             this->get_adjoint_rhs(i).close();
01045             sensitivities[i][k] += (this->get_adjoint_rhs(i).dot(this->get_weighted_sensitivity_solution()) -
01046                                     this->rhs->dot(this->get_weighted_sensitivity_adjoint_solution(i)) -
01047                                     this->get_adjoint_solution(i).dot(*tempvec)) / (2.*delta_p);
01048           }
01049 
01050       *parameters[k] = old_parameter - delta_p;
01051       this->assembly(true, true);
01052       this->rhs->close();
01053       this->matrix->close();
01054       this->assemble_qoi_derivative(qoi_indices);
01055 
01056       this->matrix->vector_mult(*tempvec, this->get_weighted_sensitivity_solution());
01057 
01058       for (unsigned int i=0; i != Nq; ++i)
01059         if (qoi_indices.has_index(i))
01060           {
01061             this->get_adjoint_rhs(i).close();
01062             sensitivities[i][k] += (-this->get_adjoint_rhs(i).dot(this->get_weighted_sensitivity_solution()) +
01063                                     this->rhs->dot(this->get_weighted_sensitivity_adjoint_solution(i)) +
01064                                     this->get_adjoint_solution(i).dot(*tempvec)) / (2.*delta_p);
01065           }
01066     }
01067 
01068   // All parameters have been reset.
01069   // Don't leave the qoi or system changed - principle of least
01070   // surprise.
01071   this->assembly(true, true);
01072   this->rhs->close();
01073   this->matrix->close();
01074   this->assemble_qoi(qoi_indices);
01075 }

void libMesh::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 509 of file system.C.

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

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

void libMesh::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.

Reimplemented in libMesh::TransientSystem< RBConstruction >.

Definition at line 431 of file system.C.

References libMesh::System::current_local_solution, libMesh::System::get_dof_map(), libMesh::DofMap::get_send_list(), libMesh::System::n_vars(), and libMesh::System::solution.

00432 {
00433   parallel_only();
00434 
00435   // If this system is empty... don't do anything!
00436   if(!this->n_vars())
00437     return;
00438 
00439   const std::vector<dof_id_type>& send_list = this->get_dof_map().get_send_list ();
00440 
00441   // Check sizes
00442   libmesh_assert_equal_to (current_local_solution->size(), solution->size());
00443   // Not true with ghosted vectors
00444   // libmesh_assert_equal_to (current_local_solution->local_size(), solution->size());
00445   // libmesh_assert (!send_list.empty());
00446   libmesh_assert_less_equal (send_list.size(), solution->size());
00447 
00448   // Create current_local_solution from solution.  This will
00449   // put a local copy of solution into current_local_solution.
00450   solution->localize (*current_local_solution, send_list);
00451 }

void libMesh::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 113 of file system_io.C.

References libMesh::System::_additional_data_written, libMesh::System::_written_var_indices, libMesh::System::add_variable(), libMesh::System::add_vector(), libMesh::Parallel::Communicator::broadcast(), libMesh::System::clear(), libMesh::CommWorld, libMesh::Xdr::data(), libMesh::FEType::family, libMesh::System::get_mesh(), libMesh::FEType::inf_map, libMesh::MeshBase::mesh_dimension(), libMeshEnums::MONOMIAL, libMesh::on_command_line(), libMesh::FEType::order, libMesh::out, libMesh::processor_id(), libMesh::FEType::radial_family, libMesh::FEType::radial_order, libMesh::Xdr::reading(), libMesh::System::variable_number(), libMesh::Xdr::version(), and libMeshEnums::XYZ.

Referenced by libMesh::EquationSystems::_read_impl().

00118 {
00119   // This method implements the input of a
00120   // System object, embedded in the output of
00121   // an EquationSystems<T_sys>.  This warrants some
00122   // documentation.  The output file essentially
00123   // consists of 5 sections:
00124   //
00125   // for this system
00126   //
00127   //   5.) The number of variables in the system (unsigned int)
00128   //
00129   //   for each variable in the system
00130   //
00131   //     6.) The name of the variable (string)
00132   //
00133   //     6.1.) Variable subdmains
00134   //
00135   //     7.) Combined in an FEType:
00136   //         - The approximation order(s) of the variable
00137   //           (Order Enum, cast to int/s)
00138   //         - The finite element family/ies of the variable
00139   //           (FEFamily Enum, cast to int/s)
00140   //
00141   //   end variable loop
00142   //
00143   //   8.) The number of additional vectors (unsigned int),
00144   //
00145   //     for each additional vector in the system object
00146   //
00147   //     9.) the name of the additional vector  (string)
00148   //
00149   // end system
00150   libmesh_assert (io.reading());
00151 
00152   // Possibly clear data structures and start from scratch.
00153   if (read_header_in)
00154     this->clear ();
00155 
00156   // Figure out if we need to read infinite element information.
00157   // This will be true if the version string contains " with infinite elements"
00158   const bool read_ifem_info =
00159     (version.rfind(" with infinite elements") < version.size()) ||
00160     libMesh::on_command_line ("--read_ifem_systems");
00161 
00162 
00163   {
00164     // 5.)
00165     // Read the number of variables in the system
00166     unsigned int nv=0;
00167     if (libMesh::processor_id() == 0)
00168       io.data (nv);
00169     CommWorld.broadcast(nv);
00170 
00171     _written_var_indices.clear();
00172     _written_var_indices.resize(nv, 0);
00173 
00174     for (unsigned int var=0; var<nv; var++)
00175       {
00176         // 6.)
00177         // Read the name of the var-th variable
00178         std::string var_name;
00179         if (libMesh::processor_id() == 0)
00180           io.data (var_name);
00181         CommWorld.broadcast(var_name);
00182 
00183         // 6.1.)
00184         std::set<subdomain_id_type> domains;
00185         if (io.version() >= LIBMESH_VERSION_ID(0,7,2))
00186         {
00187           std::vector<subdomain_id_type> domain_array;
00188           if (libMesh::processor_id() == 0)
00189             io.data (domain_array);
00190           for (std::vector<subdomain_id_type>::iterator it = domain_array.begin(); it != domain_array.end(); ++it)
00191                   domains.insert(*it);
00192         }
00193         CommWorld.broadcast(domains);
00194 
00195         // 7.)
00196         // Read the approximation order(s) of the var-th variable
00197         int order=0;
00198         if (libMesh::processor_id() == 0)
00199           io.data (order);
00200         CommWorld.broadcast(order);
00201 
00202 
00203         // do the same for infinite element radial_order
00204         int rad_order=0;
00205         if (read_ifem_info)
00206           {
00207             if (libMesh::processor_id() == 0)
00208               io.data(rad_order);
00209             CommWorld.broadcast(rad_order);
00210           }
00211 
00212 
00213         // Read the finite element type of the var-th variable
00214         int fam=0;
00215         if (libMesh::processor_id() == 0)
00216           io.data (fam);
00217         CommWorld.broadcast(fam);
00218         FEType type;
00219         type.order  = static_cast<Order>(order);
00220         type.family = static_cast<FEFamily>(fam);
00221 
00222         // Check for incompatibilities.  The shape function indexing was
00223         // changed for the monomial and xyz finite element families to
00224         // simplify extension to arbitrary p.  The consequence is that
00225         // old restart files will not be read correctly.  This is expected
00226         // to be an unlikely occurance, but catch it anyway.
00227         if (read_legacy_format)
00228           if ((type.family == MONOMIAL || type.family == XYZ) &&
00229               ((type.order > 2 && this->get_mesh().mesh_dimension() == 2) ||
00230                (type.order > 1 && this->get_mesh().mesh_dimension() == 3)))
00231             {
00232               libmesh_here();
00233               libMesh::out << "*****************************************************************\n"
00234                             << "* WARNING: reading a potentially incompatible restart file!!!   *\n"
00235                             << "*  contact libmesh-users@lists.sourceforge.net for more details *\n"
00236                             << "*****************************************************************"
00237                             << std::endl;
00238             }
00239 
00240         // Read additional information for infinite elements
00241         int radial_fam=0;
00242         int i_map=0;
00243         if (read_ifem_info)
00244           {
00245             if (libMesh::processor_id() == 0)
00246               io.data (radial_fam);
00247             CommWorld.broadcast(radial_fam);
00248             if (libMesh::processor_id() == 0)
00249               io.data (i_map);
00250             CommWorld.broadcast(i_map);
00251           }
00252 
00253 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
00254 
00255         type.radial_order  = static_cast<Order>(rad_order);
00256         type.radial_family = static_cast<FEFamily>(radial_fam);
00257         type.inf_map       = static_cast<InfMapType>(i_map);
00258 
00259 #endif
00260 
00261         if (read_header_in)
00262         {
00263           if (domains.empty())
00264             _written_var_indices[var] = this->add_variable (var_name, type);
00265           else
00266             _written_var_indices[var] = this->add_variable (var_name, type, &domains);
00267         }
00268     else
00269           _written_var_indices[var] = this->variable_number(var_name);
00270       }
00271   }
00272 
00273   // 8.)
00274   // Read the number of additional vectors.
00275   unsigned int nvecs=0;
00276   if (libMesh::processor_id() == 0)
00277     io.data (nvecs);
00278   CommWorld.broadcast(nvecs);
00279 
00280   // If nvecs > 0, this means that write_additional_data
00281   // was true when this file was written.  We will need to
00282   // make use of this fact later.
00283   if (nvecs > 0)
00284     this->_additional_data_written = true;
00285 
00286   for (unsigned int vec=0; vec<nvecs; vec++)
00287     {
00288       // 9.)
00289       // Read the name of the vec-th additional vector
00290       std::string vec_name;
00291       if (libMesh::processor_id() == 0)
00292         io.data (vec_name);
00293       CommWorld.broadcast(vec_name);
00294 
00295       if (read_additional_data)
00296         {
00297           // Systems now can handle adding post-initialization vectors
00298 //        libmesh_assert(this->_can_add_vectors);
00299           // Some systems may have added their own vectors already
00300 //        libmesh_assert_equal_to (this->_vectors.count(vec_name), 0);
00301 
00302           this->add_vector(vec_name);
00303         }
00304     }
00305 }

void libMesh::System::read_legacy_data	(	Xdr &	io,
		const bool	read_additional_data = `true`
	)			`[inherited]`

Reads additional data, namely vectors, for this System.

Definition at line 309 of file system_io.C.

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

00311 {
00312   libmesh_deprecated();
00313 
00314   // This method implements the output of the vectors
00315   // contained in this System object, embedded in the
00316   // output of an EquationSystems<T_sys>.
00317   //
00318   //   10.) The global solution vector, re-ordered to be node-major
00319   //       (More on this later.)
00320   //
00321   //      for each additional vector in the object
00322   //
00323   //      11.) The global additional vector, re-ordered to be
00324   //           node-major (More on this later.)
00325   libmesh_assert (io.reading());
00326 
00327   // read and reordering buffers
00328   std::vector<Number> global_vector;
00329   std::vector<Number> reordered_vector;
00330 
00331   // 10.)
00332   // Read and set the solution vector
00333   {
00334     if (libMesh::processor_id() == 0)
00335       io.data (global_vector);
00336     CommWorld.broadcast(global_vector);
00337 
00338     // Remember that the stored vector is node-major.
00339     // We need to put it into whatever application-specific
00340     // ordering we may have using the dof_map.
00341     reordered_vector.resize(global_vector.size());
00342 
00343     //libMesh::out << "global_vector.size()=" << global_vector.size() << std::endl;
00344     //libMesh::out << "this->n_dofs()=" << this->n_dofs() << std::endl;
00345 
00346     libmesh_assert_equal_to (global_vector.size(), this->n_dofs());
00347 
00348     dof_id_type cnt=0;
00349 
00350     const unsigned int sys = this->number();
00351     const unsigned int nv  = this->_written_var_indices.size();
00352     libmesh_assert_less_equal (nv, this->n_vars());
00353 
00354     for (unsigned int data_var=0; data_var<nv; data_var++)
00355       {
00356         const unsigned int var = _written_var_indices[data_var];
00357 
00358         // First reorder the nodal DOF values
00359         {
00360           MeshBase::node_iterator
00361             it  = this->get_mesh().nodes_begin(),
00362             end = this->get_mesh().nodes_end();
00363 
00364           for (; it != end; ++it)
00365             for (unsigned int index=0; index<(*it)->n_comp(sys,var); index++)
00366               {
00367                 libmesh_assert_not_equal_to ((*it)->dof_number(sys, var, index),
00368                                              DofObject::invalid_id);
00369 
00370                 libmesh_assert_less (cnt, global_vector.size());
00371 
00372                 reordered_vector[(*it)->dof_number(sys, var, index)] =
00373                   global_vector[cnt++];
00374               }
00375         }
00376 
00377         // Then reorder the element DOF values
00378         {
00379           MeshBase::element_iterator
00380             it  = this->get_mesh().active_elements_begin(),
00381             end = this->get_mesh().active_elements_end();
00382 
00383           for (; it != end; ++it)
00384             for (unsigned int index=0; index<(*it)->n_comp(sys,var); index++)
00385               {
00386                 libmesh_assert_not_equal_to ((*it)->dof_number(sys, var, index),
00387                                              DofObject::invalid_id);
00388 
00389                 libmesh_assert_less (cnt, global_vector.size());
00390 
00391                 reordered_vector[(*it)->dof_number(sys, var, index)] =
00392                   global_vector[cnt++];
00393               }
00394         }
00395       }
00396 
00397     *(this->solution) = reordered_vector;
00398   }
00399 
00400   // For each additional vector, simply go through the list.
00401   // ONLY attempt to do this IF additional data was actually
00402   // written to the file for this system (controlled by the
00403   // _additional_data_written flag).
00404   if (this->_additional_data_written)
00405     {
00406       std::map<std::string, NumericVector<Number>* >::iterator
00407         pos = this->_vectors.begin();
00408 
00409       for (; pos != this->_vectors.end(); ++pos)
00410         {
00411           // 11.)
00412           // Read the values of the vec-th additional vector.
00413           // Prior do _not_ clear, but fill with zero, since the
00414           // additional vectors _have_ to have the same size
00415           // as the solution vector
00416           std::fill (global_vector.begin(), global_vector.end(), libMesh::zero);
00417 
00418           if (libMesh::processor_id() == 0)
00419             io.data (global_vector);
00420           CommWorld.broadcast(global_vector);
00421 
00422           // If read_additional_data==true, then we will keep this vector, otherwise
00423           // we are going to throw it away.
00424           if (read_additional_data)
00425             {
00426               // Remember that the stored vector is node-major.
00427               // We need to put it into whatever application-specific
00428               // ordering we may have using the dof_map.
00429               std::fill (reordered_vector.begin(),
00430                          reordered_vector.end(),
00431                          libMesh::zero);
00432 
00433               reordered_vector.resize(global_vector.size());
00434 
00435               libmesh_assert_equal_to (global_vector.size(), this->n_dofs());
00436 
00437               dof_id_type cnt=0;
00438 
00439               const unsigned int sys = this->number();
00440               const unsigned int nv  = this->_written_var_indices.size();
00441               libmesh_assert_less_equal (nv, this->n_vars());
00442 
00443               for (unsigned int data_var=0; data_var<nv; data_var++)
00444                 {
00445                   const unsigned int var = _written_var_indices[data_var];
00446                   // First reorder the nodal DOF values
00447                   {
00448                     MeshBase::node_iterator
00449                       it  = this->get_mesh().nodes_begin(),
00450                       end = this->get_mesh().nodes_end();
00451 
00452                     for (; it!=end; ++it)
00453                       for (unsigned int index=0; index<(*it)->n_comp(sys,var); index++)
00454                         {
00455                           libmesh_assert_not_equal_to ((*it)->dof_number(sys, var, index),
00456                                                        DofObject::invalid_id);
00457 
00458                           libmesh_assert_less (cnt, global_vector.size());
00459 
00460                           reordered_vector[(*it)->dof_number(sys, var, index)] =
00461                             global_vector[cnt++];
00462                         }
00463                   }
00464 
00465                   // Then reorder the element DOF values
00466                   {
00467                     MeshBase::element_iterator
00468                       it  = this->get_mesh().active_elements_begin(),
00469                       end = this->get_mesh().active_elements_end();
00470 
00471                     for (; it!=end; ++it)
00472                       for (unsigned int index=0; index<(*it)->n_comp(sys,var); index++)
00473                         {
00474                           libmesh_assert_not_equal_to ((*it)->dof_number(sys, var, index),
00475                                                        DofObject::invalid_id);
00476 
00477                           libmesh_assert_less (cnt, global_vector.size());
00478 
00479                           reordered_vector[(*it)->dof_number(sys, var, index)] =
00480                             global_vector[cnt++];
00481                         }
00482                   }
00483                 }
00484 
00485               // use the overloaded operator=(std::vector) to assign the values
00486               *(pos->second) = reordered_vector;
00487             }
00488         }
00489     } // end if (_additional_data_written)
00490 }

void libMesh::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 494 of file system_io.C.

References libMesh::System::_vectors, libMesh::System::_written_var_indices, libMesh::Xdr::data(), libMesh::FEType::family, libMesh::System::get_dof_map(), libMesh::System::get_mesh(), libMesh::DofObject::invalid_id, libMesh::Xdr::is_open(), libMesh::n_processors(), libMesh::System::n_vars(), libMesh::System::number(), libMesh::processor_id(), libMesh::Xdr::reading(), libMeshEnums::SCALAR, libMesh::DofMap::SCALAR_dof_indices(), libMesh::System::solution, libMesh::Variable::type(), and libMesh::System::variable().

00496 {
00516   // PerfLog pl("IO Performance",false);
00517   // pl.push("read_parallel_data");
00518   dof_id_type total_read_size = 0;
00519   
00520   libmesh_assert (io.reading());
00521   libmesh_assert (io.is_open());
00522 
00523   // build the ordered nodes and element maps.
00524   // when writing/reading parallel files we need to iterate
00525   // over our nodes/elements in order of increasing global id().
00526   // however, this is not guaranteed to be ordering we obtain
00527   // by using the node_iterators/element_iterators directly.
00528   // so build a set, sorted by id(), that provides the ordering.
00529   // further, for memory economy build the set but then transfer
00530   // its contents to vectors, which will be sorted.
00531   std::vector<const DofObject*> ordered_nodes, ordered_elements;
00532   {
00533     std::set<const DofObject*, CompareDofObjectsByID>
00534       ordered_nodes_set (this->get_mesh().local_nodes_begin(),
00535                          this->get_mesh().local_nodes_end());
00536 
00537       ordered_nodes.insert(ordered_nodes.end(),
00538                            ordered_nodes_set.begin(),
00539                            ordered_nodes_set.end());
00540   }
00541   {
00542     std::set<const DofObject*, CompareDofObjectsByID>
00543       ordered_elements_set (this->get_mesh().local_elements_begin(),
00544                             this->get_mesh().local_elements_end());
00545 
00546       ordered_elements.insert(ordered_elements.end(),
00547                               ordered_elements_set.begin(),
00548                               ordered_elements_set.end());
00549   }
00550 
00551   std::vector<Number> io_buffer;
00552 
00553   // 9.)
00554   //
00555   // Actually read the solution components
00556   // for the ith system to disk
00557   io.data(io_buffer);
00558 
00559   total_read_size += io_buffer.size();
00560 
00561   const unsigned int sys_num = this->number();
00562   const unsigned int nv      = this->_written_var_indices.size();
00563   libmesh_assert_less_equal (nv, this->n_vars());
00564 
00565   dof_id_type cnt=0;
00566 
00567   // Loop over each non-SCALAR variable and each node, and read out the value.
00568   for (unsigned int data_var=0; data_var<nv; data_var++)
00569     {
00570       const unsigned int var = _written_var_indices[data_var];
00571       if(this->variable(var).type().family != SCALAR)
00572         {
00573           // First read the node DOF values
00574           for (std::vector<const DofObject*>::const_iterator
00575                it = ordered_nodes.begin(); it != ordered_nodes.end(); ++it)
00576             for (unsigned int comp=0; comp<(*it)->n_comp(sys_num, var); comp++)
00577               {
00578                 libmesh_assert_not_equal_to ((*it)->dof_number(sys_num, var, comp),
00579                                              DofObject::invalid_id);
00580                 libmesh_assert_less (cnt, io_buffer.size());
00581                 this->solution->set((*it)->dof_number(sys_num, var, comp), io_buffer[cnt++]);
00582               }
00583 
00584           // Then read the element DOF values
00585           for (std::vector<const DofObject*>::const_iterator
00586              it = ordered_elements.begin(); it != ordered_elements.end(); ++it)
00587             for (unsigned int comp=0; comp<(*it)->n_comp(sys_num, var); comp++)
00588             {
00589               libmesh_assert_not_equal_to ((*it)->dof_number(sys_num, var, comp),
00590                                            DofObject::invalid_id);
00591               libmesh_assert_less (cnt, io_buffer.size());
00592               this->solution->set((*it)->dof_number(sys_num, var, comp), io_buffer[cnt++]);
00593             }
00594         }
00595     }
00596 
00597   // Finally, read the SCALAR variables on the last processor
00598   for (unsigned int data_var=0; data_var<nv; data_var++)
00599     {
00600       const unsigned int var = _written_var_indices[data_var];
00601       if(this->variable(var).type().family == SCALAR)
00602         {
00603           if (libMesh::processor_id() == (libMesh::n_processors()-1))
00604             {
00605               const DofMap& dof_map = this->get_dof_map();
00606               std::vector<dof_id_type> SCALAR_dofs;
00607               dof_map.SCALAR_dof_indices(SCALAR_dofs, var);
00608 
00609               for(unsigned int i=0; i<SCALAR_dofs.size(); i++)
00610                 {
00611                   this->solution->set( SCALAR_dofs[i], io_buffer[cnt++] );
00612                 }
00613             }
00614         }
00615     }
00616 
00617   // And we're done setting solution entries
00618   this->solution->close();
00619 
00620   // Only read additional vectors if wanted
00621   if (read_additional_data)
00622     {
00623       std::map<std::string, NumericVector<Number>* >::const_iterator
00624         pos = _vectors.begin();
00625 
00626       for(; pos != this->_vectors.end(); ++pos)
00627         {
00628           cnt=0;
00629           io_buffer.clear();
00630 
00631           // 10.)
00632           //
00633           // Actually read the additional vector components
00634           // for the ith system to disk
00635           io.data(io_buffer);
00636 
00637           total_read_size += io_buffer.size();
00638 
00639           // Loop over each non-SCALAR variable and each node, and read out the value.
00640           for (unsigned int data_var=0; data_var<nv; data_var++)
00641             {
00642               const unsigned int var = _written_var_indices[data_var];
00643               if(this->variable(var).type().family != SCALAR)
00644                 {
00645                   // First read the node DOF values
00646                   for (std::vector<const DofObject*>::const_iterator
00647                        it = ordered_nodes.begin(); it != ordered_nodes.end(); ++it)
00648                     for (unsigned int comp=0; comp<(*it)->n_comp(sys_num, var); comp++)
00649                       {
00650                         libmesh_assert_not_equal_to ((*it)->dof_number(sys_num, var, comp),
00651                                                      DofObject::invalid_id);
00652                         libmesh_assert_less (cnt, io_buffer.size());
00653                         pos->second->set((*it)->dof_number(sys_num, var, comp), io_buffer[cnt++]);
00654                       }
00655 
00656                   // Then read the element DOF values
00657                   for (std::vector<const DofObject*>::const_iterator
00658                        it = ordered_elements.begin(); it != ordered_elements.end(); ++it)
00659                     for (unsigned int comp=0; comp<(*it)->n_comp(sys_num, var); comp++)
00660                       {
00661                         libmesh_assert_not_equal_to ((*it)->dof_number(sys_num, var, comp),
00662                                                      DofObject::invalid_id);
00663                         libmesh_assert_less (cnt, io_buffer.size());
00664                         pos->second->set((*it)->dof_number(sys_num, var, comp), io_buffer[cnt++]);
00665                       }
00666                 }
00667             }
00668 
00669           // Finally, read the SCALAR variables on the last processor
00670           for (unsigned int data_var=0; data_var<nv; data_var++)
00671             {
00672               const unsigned int var = _written_var_indices[data_var];
00673               if(this->variable(var).type().family == SCALAR)
00674                 {
00675                   if (libMesh::processor_id() == (libMesh::n_processors()-1))
00676                     {
00677                       const DofMap& dof_map = this->get_dof_map();
00678                       std::vector<dof_id_type> SCALAR_dofs;
00679                       dof_map.SCALAR_dof_indices(SCALAR_dofs, var);
00680 
00681                       for(unsigned int i=0; i<SCALAR_dofs.size(); i++)
00682                         {
00683                           pos->second->set( SCALAR_dofs[i], io_buffer[cnt++] );
00684                         }
00685                     }
00686                 }
00687             }
00688 
00689           // And we're done setting entries for this variable
00690           pos->second->close();
00691         }
00692     }
00693 
00694   // const Real 
00695   //   dt   = pl.get_elapsed_time(),
00696   //   rate = total_read_size*sizeof(Number)/dt; 
00697 
00698   // std::cerr << "Read " << total_read_size << " \"Number\" values\n"
00699   //        << " Elapsed time = " << dt << '\n'
00700   //        << " Rate = " << rate/1.e6 << "(MB/sec)\n\n";
00701 
00702   // pl.pop("read_parallel_data");
00703 }

void libMesh::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 707 of file system_io.C.

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

00709 {
00710   // This method implements the input of the vectors
00711   // contained in this System object, embedded in the
00712   // output of an EquationSystems<T_sys>.
00713   //
00714   //   10.) The global solution vector, re-ordered to be node-major
00715   //       (More on this later.)
00716   //
00717   //      for each additional vector in the object
00718   //
00719   //      11.) The global additional vector, re-ordered to be
00720   //          node-major (More on this later.)
00721   parallel_only();
00722   std::string comment;
00723 
00724   // PerfLog pl("IO Performance",false);
00725   // pl.push("read_serialized_data");
00726   std::size_t total_read_size = 0;
00727 
00728   // 10.)
00729   // Read the global solution vector
00730   {
00731     total_read_size +=
00732       this->read_serialized_vector(io, *this->solution);
00733 
00734     // get the comment
00735     if (libMesh::processor_id() == 0)
00736       io.comment (comment);
00737   }
00738 
00739   // 11.)
00740   // Only read additional vectors if wanted
00741   if (read_additional_data)
00742     {
00743       std::map<std::string, NumericVector<Number>* >::const_iterator
00744         pos = _vectors.begin();
00745 
00746       for(; pos != this->_vectors.end(); ++pos)
00747         {
00748           total_read_size +=
00749             this->read_serialized_vector(io, *pos->second);
00750 
00751           // get the comment
00752           if (libMesh::processor_id() == 0)
00753             io.comment (comment);
00754 
00755         }
00756     }
00757   
00758   // const Real 
00759   //   dt   = pl.get_elapsed_time(),
00760   //   rate = total_read_size*sizeof(Number)/dt; 
00761 
00762   // std::cout << "Read " << total_read_size << " \"Number\" values\n"
00763   //        << " Elapsed time = " << dt << '\n'
00764   //        << " Rate = " << rate/1.e6 << "(MB/sec)\n\n";
00765 
00766   // pl.pop("read_serialized_data");
00767 }

dof_id_type libMesh::System::read_serialized_vectors	(	Xdr &	io,
		const std::vector< NumericVector< Number > * > &	vectors
	)			const `[inherited]`

Read a number of identically distributed vectors. This method allows for optimization for the multiple vector case by only communicating the metadata once.

Definition at line 2164 of file system_io.C.

References libMesh::Xdr::data(), libMesh::FEType::family, libMesh::System::get_mesh(), libMesh::MeshBase::n_elem(), libMesh::MeshTools::n_elem(), libMesh::MeshBase::n_nodes(), n_nodes, libMesh::System::n_vars(), libMesh::processor_id(), libMesh::System::read_SCALAR_dofs(), libMesh::System::read_serialized_blocked_dof_objects(), libMesh::Xdr::reading(), libMeshEnums::SCALAR, libMesh::Variable::type(), and libMesh::System::variable().

02166 {
02167   parallel_only();
02168 
02169   // Error checking
02170 // #ifndef NDEBUG
02171 //   // In parallel we better be reading a parallel vector -- if not
02172 //   // we will not set all of its components below!!
02173 //   if (libMesh::n_processors() > 1)
02174 //     {
02175 //       libmesh_assert (vec.type() == PARALLEL ||
02176 //                    vec.type() == GHOSTED);
02177 //     }
02178 // #endif
02179 
02180   libmesh_assert (io.reading());
02181 
02182   // sizes
02183   unsigned int num_vecs=0;
02184   dof_id_type vector_length=0;
02185 
02186   if (libMesh::processor_id() == 0) 
02187     {
02188       // Get the number of vectors
02189       io.data(num_vecs);
02190       // Get the buffer size
02191       io.data(vector_length);
02192 
02193       libmesh_assert_equal_to (num_vecs, vectors.size());
02194       
02195       if (num_vecs != 0)
02196         {
02197           libmesh_assert_not_equal_to (vectors[0], 0);
02198           libmesh_assert_equal_to     (vectors[0]->size(), vector_length);
02199         }
02200     }
02201 
02202   // no need to actually communicate these.
02203   // CommWorld.broadcast(num_vecs);
02204   // CommWorld.broadcast(vector_length);
02205   
02206   // Cache these - they are not free!
02207   const dof_id_type
02208     n_nodes = this->get_mesh().n_nodes(),
02209     n_elem  = this->get_mesh().n_elem();  
02210 
02211   dof_id_type read_length = 0.;
02212   
02213   //---------------------------------
02214   // Collect the values for all nodes
02215   read_length +=
02216     this->read_serialized_blocked_dof_objects (n_nodes,
02217                                                this->get_mesh().local_nodes_begin(),
02218                                                this->get_mesh().local_nodes_end(),
02219                                                io,
02220                                                vectors);
02221 
02222   //------------------------------------
02223   // Collect the values for all elements
02224   read_length +=
02225     this->read_serialized_blocked_dof_objects (n_elem,
02226                                                this->get_mesh().local_elements_begin(),
02227                                                this->get_mesh().local_elements_end(),
02228                                                io,
02229                                                vectors);
02230 
02231   //-------------------------------------------
02232   // Finally loop over all the SCALAR variables
02233   for (unsigned int vec=0; vec<vectors.size(); vec++)
02234     for (unsigned int var=0; var<this->n_vars(); var++)
02235       if(this->variable(var).type().family == SCALAR)
02236         {
02237           libmesh_assert_not_equal_to (vectors[vec], 0);
02238           
02239           read_length +=
02240             this->read_SCALAR_dofs (var, io, *vectors[vec]);
02241         }
02242 
02243   //---------------------------------------
02244   // last step - must close all the vectors
02245   for (unsigned int vec=0; vec<vectors.size(); vec++)
02246     {
02247       libmesh_assert_not_equal_to (vectors[vec], 0);      
02248       vectors[vec]->close();
02249     }
02250       
02251   return read_length;
02252 }

void libMesh::LinearImplicitSystem::reinit ( ) [virtual]

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

Reimplemented from libMesh::ImplicitSystem.

Reimplemented in libMesh::NewmarkSystem, and libMesh::TransientSystem< RBConstruction >.

Definition at line 85 of file linear_implicit_system.C.

References linear_solver, and libMesh::ImplicitSystem::reinit().

00086 {
00087   // re-initialize the linear solver interface
00088   linear_solver->clear();
00089 
00090   // initialize parent data
00091   Parent::reinit();
00092 }

void libMesh::LinearImplicitSystem::release_linear_solver ( LinearSolver< Number > * ) const [virtual]

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

Reimplemented from libMesh::ImplicitSystem.

Definition at line 367 of file linear_implicit_system.C.

00368 {
00369 }

void libMesh::System::remove_vector ( const std::string & vec_name ) [inherited]

Removes the additional vector vec_name from this system

Definition at line 711 of file system.C.

References libMesh::System::_vector_projections, libMesh::System::_vector_types, libMesh::System::_vectors, and libMesh::System::have_vector().

00712 {
00713   //Return if the vector does not exist
00714   if ( !(this->have_vector(vec_name)) )
00715     return;
00716 
00717   _vectors[vec_name]->clear();
00718   delete _vectors[vec_name];
00719   _vectors[vec_name] = NULL;
00720 
00721   _vectors.erase(vec_name);
00722   _vector_projections.erase(vec_name);
00723   _vector_types.erase(vec_name);
00724 }

SparseMatrix< Number > * libMesh::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 245 of file implicit_system.C.

References libMesh::ImplicitSystem::_matrices.

00246 {
00247   // Make sure the matrix exists
00248   matrices_iterator pos = _matrices.find (mat_name);
00249 
00250   if (pos == _matrices.end())
00251     return NULL;
00252 
00253   return pos->second;
00254 }

const SparseMatrix< Number > * libMesh::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 232 of file implicit_system.C.

References libMesh::ImplicitSystem::_matrices.

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

00233 {
00234   // Make sure the matrix exists
00235   const_matrices_iterator pos = _matrices.find (mat_name);
00236 
00237   if (pos == _matrices.end())
00238     return NULL;
00239 
00240   return pos->second;
00241 }

NumericVector< Number > * libMesh::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 767 of file system.C.

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

00768 {
00769   vectors_iterator v = vectors_begin();
00770   vectors_iterator v_end = vectors_end();
00771   unsigned int num = 0;
00772   while((num<vec_num) && (v!=v_end))
00773     {
00774       num++;
00775       ++v;
00776     }
00777   if (v==v_end)
00778     return NULL;
00779   return v->second;
00780 }

const NumericVector< Number > * libMesh::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 750 of file system.C.

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

00751 {
00752   const_vectors_iterator v = vectors_begin();
00753   const_vectors_iterator v_end = vectors_end();
00754   unsigned int num = 0;
00755   while((num<vec_num) && (v!=v_end))
00756     {
00757       num++;
00758       ++v;
00759     }
00760   if (v==v_end)
00761     return NULL;
00762   return v->second;
00763 }

NumericVector< Number > * libMesh::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 738 of file system.C.

References libMesh::System::_vectors.

00739 {
00740   vectors_iterator pos = _vectors.find(vec_name);
00741 
00742   if (pos == _vectors.end())
00743     return NULL;
00744 
00745   return pos->second;
00746 }

const NumericVector< Number > * libMesh::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 726 of file system.C.

References libMesh::System::_vectors.

00727 {
00728   const_vectors_iterator pos = _vectors.find(vec_name);
00729 
00730   if (pos == _vectors.end())
00731     return NULL;
00732 
00733   return pos->second;
00734 }

void libMesh::LinearImplicitSystem::restrict_solve_to	(	const SystemSubset *	subset,
		const SubsetSolveMode	subset_solve_mode = `SUBSET_ZERO`
	)			`[virtual]`

After calling this method, any solve will be limited to the given subset. To disable this mode, call this method with subset being a NULL pointer.

Reimplemented from libMesh::System.

Definition at line 96 of file linear_implicit_system.C.

References _subset, _subset_solve_mode, and libMesh::SystemSubset::get_system().

Referenced by clear().

00098 {
00099   _subset = subset;
00100   _subset_solve_mode = subset_solve_mode;
00101   if(subset!=NULL)
00102     {
00103       libmesh_assert_equal_to (&subset->get_system(), this);
00104     }
00105 }

void libMesh::System::restrict_vectors ( ) [virtual, inherited]

Restrict vectors after the mesh has coarsened

Definition at line 316 of file system.C.

References libMesh::System::_dof_map, libMesh::System::_solution_projection, libMesh::System::_vector_projections, libMesh::System::_vector_types, libMesh::System::_vectors, libMesh::System::current_local_solution, libMesh::err, libMeshEnums::GHOSTED, libMesh::System::n_dofs(), libMesh::System::n_local_dofs(), libMesh::System::project_vector(), and libMesh::System::solution.

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

00317 {
00318 #ifdef LIBMESH_ENABLE_AMR
00319   // Restrict the _vectors on the coarsened cells
00320   for (vectors_iterator pos = _vectors.begin(); pos != _vectors.end(); ++pos)
00321     {
00322       NumericVector<Number>* v = pos->second;
00323 
00324       if (_vector_projections[pos->first])
00325         this->project_vector (*v);
00326       else
00327       {
00328         ParallelType type = _vector_types[pos->first];
00329 
00330         if(type == GHOSTED)
00331         {
00332 #ifdef LIBMESH_ENABLE_GHOSTED
00333           pos->second->init (this->n_dofs(), this->n_local_dofs(),
00334                          _dof_map->get_send_list(), false,
00335                          GHOSTED);
00336 #else
00337           libMesh::err << "Cannot initialize ghosted vectors when they are not enabled." << std::endl;
00338           libmesh_error();
00339 #endif
00340         }
00341         else
00342           pos->second->init (this->n_dofs(), this->n_local_dofs(), false, type);
00343       }
00344     }
00345 
00346   const std::vector<dof_id_type>& send_list = _dof_map->get_send_list ();
00347 
00348   // Restrict the solution on the coarsened cells
00349   if (_solution_projection)
00350     this->project_vector (*solution);
00351 
00352 #ifdef LIBMESH_ENABLE_GHOSTED
00353   current_local_solution->init(this->n_dofs(),
00354                                this->n_local_dofs(), send_list,
00355                                false, GHOSTED);
00356 #else
00357   current_local_solution->init(this->n_dofs());
00358 #endif
00359 
00360   if (_solution_projection)
00361     solution->localize (*current_local_solution, send_list);
00362 
00363 #endif // LIBMESH_ENABLE_AMR
00364 }

std::pair< unsigned int, Real > libMesh::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 libMesh::System.

Definition at line 321 of file implicit_system.C.

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

00322 {
00323   // Log how long the linear solve takes.
00324   START_LOG("sensitivity_solve()", "ImplicitSystem");
00325 
00326   // The forward system should now already be solved.
00327   // Now assemble the corresponding sensitivity system.
00328 
00329   if (this->assemble_before_solve)
00330     {
00331       // Build the Jacobian
00332       this->assembly(false, true);
00333       this->matrix->close();
00334 
00335       // Reset and build the RHS from the residual derivatives
00336       this->assemble_residual_derivatives(parameters);
00337     }
00338 
00339   // The sensitivity problem is linear
00340   LinearSolver<Number> *linear_solver = this->get_linear_solver();
00341 
00342   // Our iteration counts and residuals will be sums of the individual
00343   // results
00344   std::pair<unsigned int, Real> solver_params =
00345     this->get_linear_solve_parameters();
00346   std::pair<unsigned int, Real> totalrval = std::make_pair(0,0.0);
00347 
00348   // Solve the linear system.
00349   SparseMatrix<Number> *pc = this->request_matrix("Preconditioner");
00350   for (unsigned int p=0; p != parameters.size(); ++p)
00351     {
00352       std::pair<unsigned int, Real> rval =
00353         linear_solver->solve (*matrix, pc,
00354                               this->add_sensitivity_solution(p),
00355                               this->get_sensitivity_rhs(p),
00356                               solver_params.second,
00357                               solver_params.first);
00358 
00359       totalrval.first  += rval.first;
00360       totalrval.second += rval.second;
00361     }
00362 
00363   // The linear solver may not have fit our constraints exactly
00364 #ifdef LIBMESH_ENABLE_CONSTRAINTS
00365   for (unsigned int p=0; p != parameters.size(); ++p)
00366     this->get_dof_map().enforce_constraints_exactly
00367       (*this, &this->get_sensitivity_solution(p),
00368        /* homogeneous = */ true);
00369 #endif
00370 
00371   this->release_linear_solver(linear_solver);
00372 
00373   // Stop logging the nonlinear solve
00374   STOP_LOG("sensitivity_solve()", "ImplicitSystem");
00375 
00376   return totalrval;
00377 }

void libMesh::System::set_adjoint_already_solved ( bool setting ) [inline, inherited]

Setter for the adjoint_already_solved boolean

Definition at line 365 of file system.h.

References libMesh::System::adjoint_already_solved.

00366   { adjoint_already_solved = setting;}

void libMesh::System::set_basic_system_only ( ) [inline, inherited]

Sets the system to be "basic only": i.e. advanced system components such as ImplicitSystem matrices may not be initialized. This is useful for efficiency in certain utility programs that never use System::solve(). This method must be called after the System or derived class is created but before it is initialized; e.g. from within EquationSystems::read()

Definition at line 1890 of file system.h.

References libMesh::System::_basic_system_only.

Referenced by libMesh::EquationSystems::_read_impl().

01891 {
01892   _basic_system_only = true;
01893 }

void libMesh::System::set_vector_preservation	(	const std::string &	vec_name,
		bool	preserve
	)			`[inherited]`

Allows one to set the boolean controlling whether the vector identified by vec_name should be "preserved": projected to new meshes, saved, etc.

Definition at line 887 of file system.C.

References libMesh::System::_vector_projections.

00889 {
00890   _vector_projections[vec_name] = preserve;
00891 }

void libMesh::LinearImplicitSystem::solve ( ) [virtual]

Assembles & solves the linear system A*x=b.

Reimplemented from libMesh::ExplicitSystem.

Reimplemented in libMesh::FrequencySystem.

Definition at line 109 of file linear_implicit_system.C.

References _final_linear_residual, _n_linear_iterations, _shell_matrix, _subset, _subset_solve_mode, assemble(), libMesh::System::assemble_before_solve, libMesh::SystemSubset::dof_ids(), libMesh::Parameters::get(), libMesh::System::get_equation_systems(), linear_solver, libMesh::ImplicitSystem::matrix, libMesh::EquationSystems::parameters, libMesh::Real, libMesh::ImplicitSystem::request_matrix(), libMesh::ExplicitSystem::rhs, libMesh::System::solution, and libMesh::System::update().

00110 {
00111   if (this->assemble_before_solve)
00112     // Assemble the linear system
00113     this->assemble ();
00114 
00115   // Log how long the linear solve takes.
00116   // This gets done by the LinearSolver classes now [RHS]
00117   // START_LOG("solve()", "System");
00118 
00119   // Get a reference to the EquationSystems
00120   const EquationSystems& es =
00121     this->get_equation_systems();
00122 
00123   // Get the user-specifiied linear solver tolerance
00124   const Real tol            =
00125     es.parameters.get<Real>("linear solver tolerance");
00126 
00127   // Get the user-specified maximum # of linear solver iterations
00128   const unsigned int maxits =
00129     es.parameters.get<unsigned int>("linear solver maximum iterations");
00130 
00131   if(_subset!=NULL)
00132     {
00133       linear_solver->restrict_solve_to(&_subset->dof_ids(),_subset_solve_mode);
00134     }
00135 
00136   // Solve the linear system.  Several cases:
00137   std::pair<unsigned int, Real> rval = std::make_pair(0,0.0);
00138   if(_shell_matrix)
00139     // 1.) Shell matrix with or without user-supplied preconditioner.
00140     rval = linear_solver->solve(*_shell_matrix, this->request_matrix("Preconditioner"), *solution, *rhs, tol, maxits);
00141   else
00142     // 2.) No shell matrix, with or without user-supplied preconditioner
00143     rval = linear_solver->solve (*matrix, this->request_matrix("Preconditioner"), *solution, *rhs, tol, maxits);
00144 
00145   if(_subset!=NULL)
00146     {
00147       linear_solver->restrict_solve_to(NULL);
00148     }
00149 
00150   // Store the number of linear iterations required to
00151   // solve and the final residual.
00152   _n_linear_iterations   = rval.first;
00153   _final_linear_residual = rval.second;
00154 
00155   // Stop logging the linear solve
00156   // This gets done by the LinearSolver classes now [RHS]
00157   // STOP_LOG("solve()", "System");
00158 
00159   // Update the system after the solve
00160   this->update();
00161 }

sys_type& libMesh::LinearImplicitSystem::system ( ) [inline]

Returns:: a clever pointer to the system.

Reimplemented from libMesh::ImplicitSystem.

Reimplemented in libMesh::RBConstruction, libMesh::RBConstructionBase< LinearImplicitSystem >, and libMesh::TransientSystem< RBConstruction >.

Definition at line 79 of file linear_implicit_system.h.

00079 { return *this; }

virtual std::string libMesh::LinearImplicitSystem::system_type ( ) const [inline, virtual]

Returns:: "LinearImplicit". Helps in identifying the system type in an equation system file.

Reimplemented from libMesh::ImplicitSystem.

Reimplemented in libMesh::FrequencySystem, libMesh::NewmarkSystem, libMesh::RBConstruction, and libMesh::TransientSystem< RBConstruction >.

Definition at line 141 of file linear_implicit_system.h.

00141 { return "LinearImplicit"; }

void libMesh::System::update ( ) [virtual, inherited]

Update the local values to reflect the solution on neighboring processors.

Definition at line 410 of file system.C.

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

Referenced by libMesh::__libmesh_petsc_diff_solver_jacobian(), libMesh::__libmesh_petsc_diff_solver_residual(), libMesh::FEMSystem::assemble_qoi(), libMesh::FEMSystem::assemble_qoi_derivative(), libMesh::NonlinearImplicitSystem::assembly(), libMesh::FEMSystem::assembly(), libMesh::Problem_Interface::computeF(), libMesh::Problem_Interface::computeJacobian(), libMesh::Problem_Interface::computePreconditioner(), libMesh::GMVIO::copy_nodal_solution(), DMFunction_libMesh(), DMJacobian_libMesh(), libMesh::AdjointRefinementEstimator::estimate_error(), libMesh::CondensedEigenSystem::get_eigenpair(), libMesh::FEMSystem::mesh_position_get(), libMesh::FEMSystem::postprocess(), libMesh::NonlinearImplicitSystem::solve(), libMesh::NewtonSolver::solve(), solve(), and libMesh::ExplicitSystem::solve().

00411 {
00412   libmesh_assert(solution->closed());
00413 
00414   const std::vector<dof_id_type>& send_list = _dof_map->get_send_list ();
00415 
00416   // Check sizes
00417   libmesh_assert_equal_to (current_local_solution->size(), solution->size());
00418 // More processors than elements => empty send_list
00419 //  libmesh_assert (!send_list.empty());
00420   libmesh_assert_less_equal (send_list.size(), solution->size());
00421 
00422   // Create current_local_solution from solution.  This will
00423   // put a local copy of solution into current_local_solution.
00424   // Only the necessary values (specified by the send_list)
00425   // are copied to minimize communication
00426   solution->localize (*current_local_solution, send_list);
00427 }

void libMesh::System::update_global_solution	(	std::vector< Number > &	global_soln,
		const unsigned int	dest_proc
	)			const `[inherited]`

Fill the input vector global_soln so that it contains the global solution on processor dest_proc. Requires communication with all other processors.

Definition at line 665 of file system.C.

References libMesh::System::solution.

00667 {
00668   global_soln.resize        (solution->size());
00669 
00670   solution->localize_to_one (global_soln, dest_proc);
00671 }

void libMesh::System::update_global_solution ( std::vector< Number > & global_soln ) const [inherited]

Fill the input vector global_soln so that it contains the global solution on all processors. Requires communication with all other processors.

Definition at line 656 of file system.C.

References libMesh::System::solution.

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

00657 {
00658   global_soln.resize (solution->size());
00659 
00660   solution->localize (global_soln);
00661 }

void libMesh::System::user_assembly ( ) [virtual, inherited]

Calls user's attached assembly function, or is overloaded by the user in derived classes.

Definition at line 1909 of file system.C.

References libMesh::System::_assemble_system_function, libMesh::System::_assemble_system_object, libMesh::System::_equation_systems, libMesh::System::Assembly::assemble(), and libMesh::System::name().

Referenced by libMesh::System::assemble().

01910 {
01911   // Call the user-provided assembly function,
01912   // if it was provided
01913   if (_assemble_system_function != NULL)
01914     this->_assemble_system_function (_equation_systems, this->name());
01915 
01916   // ...or the user-provided assembly object.
01917   else if (_assemble_system_object != NULL)
01918     this->_assemble_system_object->assemble();
01919 }

void libMesh::System::user_constrain ( ) [virtual, inherited]

Calls user's attached constraint function, or is overloaded by the user in derived classes.

Definition at line 1923 of file system.C.

References libMesh::System::_constrain_system_function, libMesh::System::_constrain_system_object, libMesh::System::_equation_systems, libMesh::System::Constraint::constrain(), and libMesh::System::name().

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

01924 {
01925   // Call the user-provided constraint function,
01926   // if it was provided
01927   if (_constrain_system_function!= NULL)
01928     this->_constrain_system_function(_equation_systems, this->name());
01929 
01930   // ...or the user-provided constraint object.
01931   else if (_constrain_system_object != NULL)
01932     this->_constrain_system_object->constrain();
01933 }

void libMesh::System::user_initialization ( ) [virtual, inherited]

Calls user's attached initialization function, or is overloaded by the user in derived classes.

Definition at line 1895 of file system.C.

References libMesh::System::_equation_systems, libMesh::System::_init_system_function, libMesh::System::_init_system_object, libMesh::System::Initialization::initialize(), and libMesh::System::name().

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

01896 {
01897   // Call the user-provided intialization function,
01898   // if it was provided
01899   if (_init_system_function != NULL)
01900     this->_init_system_function (_equation_systems, this->name());
01901 
01902   // ...or the user-provided initialization object.
01903   else if (_init_system_object != NULL)
01904     this->_init_system_object->initialize();
01905 }

void libMesh::System::user_QOI ( const QoISet & qoi_indices ) [virtual, inherited]

Calls user's attached quantity of interest function, or is overloaded by the user in derived classes.

Definition at line 1937 of file system.C.

References libMesh::System::_equation_systems, libMesh::System::_qoi_evaluate_function, libMesh::System::_qoi_evaluate_object, libMesh::System::name(), and libMesh::System::QOI::qoi().

Referenced by libMesh::System::assemble_qoi().

01938 {
01939   // Call the user-provided quantity of interest function,
01940   // if it was provided
01941   if (_qoi_evaluate_function != NULL)
01942     this->_qoi_evaluate_function(_equation_systems, this->name(), qoi_indices);
01943 
01944   // ...or the user-provided QOI function object.
01945   else if (_qoi_evaluate_object != NULL)
01946     this->_qoi_evaluate_object->qoi(qoi_indices);
01947 }

void libMesh::System::user_QOI_derivative ( const QoISet & qoi_indices ) [virtual, inherited]

Calls user's attached quantity of interest derivative function, or is overloaded by the user in derived classes.

Definition at line 1951 of file system.C.

References libMesh::System::_equation_systems, libMesh::System::_qoi_evaluate_derivative_function, libMesh::System::_qoi_evaluate_derivative_object, libMesh::System::name(), and libMesh::System::QOIDerivative::qoi_derivative().

Referenced by libMesh::System::assemble_qoi_derivative().

01952 {
01953   // Call the user-provided quantity of interest derivative,
01954   // if it was provided
01955   if (_qoi_evaluate_derivative_function != NULL)
01956     this->_qoi_evaluate_derivative_function(_equation_systems, this->name(), qoi_indices);
01957 
01958   // ...or the user-provided QOI derivative function object.
01959   else if (_qoi_evaluate_derivative_object != NULL)
01960     this->_qoi_evaluate_derivative_object->qoi_derivative(qoi_indices);
01961 }

const Variable & libMesh::System::variable ( unsigned int var ) const [inline, inherited]

Return a constant reference to Variable var.

Definition at line 1926 of file system.h.

References libMesh::System::_variables.

Referenced by libMesh::EquationSystems::build_solution_vector(), libMesh::WrappedFunction< Output >::component(), libMesh::EquationSystems::get_solution(), libMesh::WrappedFunction< Output >::operator()(), libMesh::System::project_vector(), libMesh::System::read_parallel_data(), libMesh::System::read_SCALAR_dofs(), libMesh::System::read_serialized_vector(), libMesh::System::read_serialized_vectors(), libMesh::System::write_header(), libMesh::System::write_parallel_data(), libMesh::System::write_serialized_vector(), and libMesh::System::write_serialized_vectors().

01927 {
01928   libmesh_assert_less (i, _variables.size());
01929 
01930   return _variables[i];
01931 }

const VariableGroup & libMesh::System::variable_group ( unsigned int vg ) const [inline, inherited]

Return a constant reference to VariableGroup vg.

Definition at line 1936 of file system.h.

References libMesh::System::_variable_groups.

Referenced by libMesh::System::get_info().

01937 {
01938   libmesh_assert_less (vg, _variable_groups.size());
01939 
01940   return _variable_groups[vg];
01941 }

const std::string & libMesh::System::variable_name ( const unsigned int i ) const [inline, inherited]

Returns:: the name of variable i.

Definition at line 1946 of file system.h.

References libMesh::System::_variables.

Referenced by libMesh::System::add_variable(), libMesh::System::add_variables(), libMesh::KellyErrorEstimator::boundary_side_integration(), libMesh::DiscontinuityMeasure::boundary_side_integration(), libMesh::WrappedFunction< Output >::component(), libMesh::ExactErrorEstimator::estimate_error(), libMesh::ExactSolution::ExactSolution(), libMesh::EquationSystems::get_solution(), libMesh::WrappedFunction< Output >::operator()(), and libMesh::System::write_header().

01947 {
01948   libmesh_assert_less (i, _variables.size());
01949 
01950   return _variables[i].name();
01951 }

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

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

Definition at line 1239 of file system.C.

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

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

01240 {
01241   // Make sure the variable exists
01242   std::map<std::string, unsigned short int>::const_iterator
01243     pos = _variable_numbers.find(var);
01244 
01245   if (pos == _variable_numbers.end())
01246     {
01247       libMesh::err << "ERROR: variable "
01248                     << var
01249                     << " does not exist in this system!"
01250                     << std::endl;
01251       libmesh_error();
01252     }
01253   libmesh_assert_equal_to (_variables[pos->second].name(), var);
01254 
01255   return pos->second;
01256 }

unsigned int libMesh::System::variable_scalar_number	(	unsigned int	var_num,
		unsigned int	component
	)			const `[inline, inherited]`

Returns:: an index, starting from 0 for the first component of the first variable, and incrementing for each component of each (potentially vector-valued) variable in the system in order. For systems with only scalar-valued variables, this will be the same as var_num

Irony: currently our only non-scalar-valued variable type is SCALAR.

Definition at line 1967 of file system.h.

References libMesh::System::_variables.

01969 {
01970   return _variables[var_num].first_scalar_number() + component;
01971 }

unsigned int libMesh::System::variable_scalar_number	(	const std::string &	var,
		unsigned int	component
	)			const `[inline, inherited]`

Returns:: an index, starting from 0 for the first component of the first variable, and incrementing for each component of each (potentially vector-valued) variable in the system in order. For systems with only scalar-valued variables, this will be the same as variable_number(var)

Irony: currently our only non-scalar-valued variable type is SCALAR.

Definition at line 1957 of file system.h.

References libMesh::System::variable_number().

Referenced by libMesh::ExactSolution::_compute_error(), libMesh::WrappedFunction< Output >::component(), libMesh::ExactErrorEstimator::find_squared_element_error(), libMesh::WrappedFunction< Output >::operator()(), libMesh::BoundaryProjectSolution::operator()(), libMesh::ProjectFEMSolution::operator()(), libMesh::ProjectSolution::operator()(), and libMesh::System::project_vector().

01959 {
01960   return variable_scalar_number(this->variable_number(var), component);
01961 }

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

Returns:: the finite element type for variable var.

Definition at line 1986 of file system.h.

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

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

const FEType & libMesh::System::variable_type ( const unsigned int i ) const [inline, inherited]

Returns:: the finite element type variable number i.

Definition at line 1976 of file system.h.

References libMesh::System::_variables.

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

01977 {
01978   libmesh_assert_less (i, _variables.size());
01979 
01980   return _variables[i].type();
01981 }

const std::string & libMesh::System::vector_name ( const NumericVector< Number > & vec_reference ) const [inherited]

Returns:: the name of a system vector, given a reference to that vector

Definition at line 868 of file system.C.

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

00869 {
00870   const_vectors_iterator v = vectors_begin();
00871   const_vectors_iterator v_end = vectors_end();
00872 
00873   for(; v != v_end; ++v)
00874     {
00875       // Check if the current vector is the one whose name we want
00876       if(&vec_reference == v->second)
00877         break; // exit loop if it is
00878     }
00879 
00880   // Before returning, make sure we didnt loop till the end and not find any match
00881   libmesh_assert (v != v_end);
00882 
00883   // Return the string associated with the current vector
00884   return v->first;
00885 }

const std::string & libMesh::System::vector_name ( const unsigned int vec_num ) const [inherited]

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

Definition at line 854 of file system.C.

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

00855 {
00856   const_vectors_iterator v = vectors_begin();
00857   const_vectors_iterator v_end = vectors_end();
00858   unsigned int num = 0;
00859   while((num<vec_num) && (v!=v_end))
00860     {
00861       num++;
00862       ++v;
00863     }
00864   libmesh_assert (v != v_end);
00865   return v->first;
00866 }

bool libMesh::System::vector_preservation ( const std::string & vec_name ) const [inherited]

Returns:: the boolean describing whether the vector identified by vec_name should be "preserved": projected to new meshes, saved, etc.

Definition at line 895 of file system.C.

References libMesh::System::_vector_projections.

Referenced by libMesh::MemorySolutionHistory::store().

00896 {
00897   if (_vector_projections.find(vec_name) == _vector_projections.end())
00898     return false;
00899 
00900   return _vector_projections.find(vec_name)->second;
00901 }

System::const_vectors_iterator libMesh::System::vectors_begin ( ) const [inline, inherited]

Beginning of vectors container

Definition at line 2044 of file system.h.

References libMesh::System::_vectors.

02045 {
02046   return _vectors.begin();
02047 }

System::vectors_iterator libMesh::System::vectors_begin ( ) [inline, inherited]

Beginning of vectors container

Definition at line 2038 of file system.h.

References libMesh::System::_vectors.

Referenced by libMesh::AdjointRefinementEstimator::estimate_error(), libMesh::System::get_vector(), libMesh::System::request_vector(), libMesh::MemorySolutionHistory::store(), libMesh::VTKIO::system_vectors_to_vtk(), and libMesh::System::vector_name().

02039 {
02040   return _vectors.begin();
02041 }

System::const_vectors_iterator libMesh::System::vectors_end ( ) const [inline, inherited]

End of vectors container

Definition at line 2056 of file system.h.

References libMesh::System::_vectors.

02057 {
02058   return _vectors.end();
02059 }

System::vectors_iterator libMesh::System::vectors_end ( ) [inline, inherited]

End of vectors container

Definition at line 2050 of file system.h.

References libMesh::System::_vectors.

Referenced by libMesh::AdjointRefinementEstimator::estimate_error(), libMesh::System::get_vector(), libMesh::System::request_vector(), libMesh::MemorySolutionHistory::store(), libMesh::VTKIO::system_vectors_to_vtk(), and libMesh::System::vector_name().

02051 {
02052   return _vectors.end();
02053 }

std::pair< unsigned int, Real > libMesh::ImplicitSystem::weighted_sensitivity_adjoint_solve	(	const ParameterVector &	parameters,
		const ParameterVector &	weights,
		const QoISet &	qoi_indices = `QoISet()`
	)			`[virtual, inherited]`

Assembles & solves the linear system(s) (dR/du)^T*z_w = sum(w_p*(d^2q/dudp - d^2R/dudp*z)), for those parameters p contained within parameters, weighted by the values w_p found within weights.

Assumes that adjoint_solve has already calculated z for each qoi in qoi_indices.

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

Reimplemented from libMesh::System.

Definition at line 437 of file implicit_system.C.

References libMesh::System::add_weighted_sensitivity_adjoint_solution(), libMesh::ExplicitSystem::assemble_qoi_derivative(), libMesh::ImplicitSystem::assembly(), libMesh::NumericVector< T >::close(), libMesh::SparseMatrix< T >::close(), libMesh::ParameterVector::deep_copy(), libMesh::DofMap::enforce_constraints_exactly(), libMesh::System::get_adjoint_rhs(), libMesh::System::get_adjoint_solution(), libMesh::System::get_dof_map(), libMesh::ImplicitSystem::get_linear_solve_parameters(), libMesh::ImplicitSystem::get_linear_solver(), libMesh::SparseMatrix< T >::get_transpose(), libMesh::System::get_weighted_sensitivity_adjoint_solution(), libMesh::QoISet::has_index(), libMesh::ImplicitSystem::matrix, libMesh::System::qoi, libMesh::Real, libMesh::AutoPtr< Tp >::release(), libMesh::ImplicitSystem::release_linear_solver(), libMesh::ExplicitSystem::rhs, libMesh::LinearSolver< T >::solve(), libMesh::TOLERANCE, libMesh::ParameterVector::value_copy(), libMesh::SparseMatrix< T >::vector_mult_add(), and libMesh::NumericVector< T >::zero_clone().

00440 {
00441   // Log how long the linear solve takes.
00442   START_LOG("weighted_sensitivity_adjoint_solve()", "ImplicitSystem");
00443 
00444   // We currently get partial derivatives via central differencing
00445   const Real delta_p = TOLERANCE;
00446 
00447   // The forward system should now already be solved.
00448   // The adjoint system should now already be solved.
00449   // Now we're assembling a weighted sum of adjoint-adjoint systems:
00450   //
00451   // dR/du (u, sum_l(w_l*z^l)) = sum_l(w_l*(Q''_ul - R''_ul (u, z)))
00452 
00453   // We'll assemble the rhs first, because the R'' term will require
00454   // perturbing the jacobian
00455 
00456   // We'll use temporary rhs vectors, because we haven't (yet) found
00457   // any good reasons why users might want to save these:
00458 
00459   std::vector<NumericVector<Number> *> temprhs(this->qoi.size());
00460   for (unsigned int i=0; i != this->qoi.size(); ++i)
00461     if (qoi_indices.has_index(i))
00462       temprhs[i] = this->rhs->zero_clone().release();
00463 
00464   // We approximate the _l partial derivatives via a central
00465   // differencing perturbation in the w_l direction:
00466   //
00467   // sum_l(w_l*v_l) ~= (v(p + dp*w_l*e_l) - v(p - dp*w_l*e_l))/(2*dp)
00468 
00469   // PETSc doesn't implement SGEMX, so neither does NumericVector,
00470   // so we want to avoid calculating f -= R'*z.  We'll thus evaluate
00471   // the above equation by first adding -v(p+dp...), then multiplying
00472   // the intermediate result vectors by -1, then adding -v(p-dp...),
00473   // then finally dividing by 2*dp.
00474 
00475   ParameterVector oldparameters, parameterperturbation;
00476   parameters.deep_copy(oldparameters);
00477   weights.deep_copy(parameterperturbation);
00478   parameterperturbation *= delta_p;
00479   parameters += parameterperturbation;
00480 
00481   this->assembly(false, true);
00482   this->matrix->close();
00483 
00484   // Take the discrete adjoint, so that we can calculate R_u(u,z) with
00485   // a matrix-vector product of R_u and z.
00486   matrix->get_transpose(*matrix);
00487 
00488   this->assemble_qoi_derivative(qoi_indices);
00489   for (unsigned int i=0; i != this->qoi.size(); ++i)
00490     if (qoi_indices.has_index(i))
00491       {
00492         this->get_adjoint_rhs(i).close();
00493         *(temprhs[i]) -= this->get_adjoint_rhs(i);
00494         this->matrix->vector_mult_add(*(temprhs[i]), this->get_adjoint_solution(i));
00495         *(temprhs[i]) *= -1.0;
00496       }
00497 
00498   oldparameters.value_copy(parameters);
00499   parameterperturbation *= -1.0;
00500   parameters += parameterperturbation;
00501 
00502   this->assembly(false, true);
00503   this->matrix->close();
00504   matrix->get_transpose(*matrix);
00505 
00506   this->assemble_qoi_derivative(qoi_indices);
00507   for (unsigned int i=0; i != this->qoi.size(); ++i)
00508     if (qoi_indices.has_index(i))
00509       {
00510         this->get_adjoint_rhs(i).close();
00511         *(temprhs[i]) -= this->get_adjoint_rhs(i);
00512         this->matrix->vector_mult_add(*(temprhs[i]), this->get_adjoint_solution(i));
00513         *(temprhs[i]) /= (2.0*delta_p);
00514       }
00515 
00516   // Finally, assemble the jacobian at the non-perturbed parameter
00517   // values.  Ignore assemble_before_solve; if we had a good
00518   // non-perturbed matrix before we've already overwritten it.
00519   oldparameters.value_copy(parameters);
00520 
00521   // if (this->assemble_before_solve)
00522     {
00523       // Build the Jacobian
00524       this->assembly(false, true);
00525       this->matrix->close();
00526 
00527       // Take the discrete adjoint
00528       matrix->get_transpose(*matrix);
00529     }
00530 
00531   // The weighted adjoint-adjoint problem is linear
00532   LinearSolver<Number> *linear_solver = this->get_linear_solver();
00533 
00534   // Our iteration counts and residuals will be sums of the individual
00535   // results
00536   std::pair<unsigned int, Real> solver_params =
00537     this->get_linear_solve_parameters();
00538   std::pair<unsigned int, Real> totalrval = std::make_pair(0,0.0);
00539 
00540   for (unsigned int i=0; i != this->qoi.size(); ++i)
00541     if (qoi_indices.has_index(i))
00542       {
00543         const std::pair<unsigned int, Real> rval =
00544           linear_solver->solve (*matrix, this->add_weighted_sensitivity_adjoint_solution(i),
00545                                  *(temprhs[i]),
00546                                  solver_params.second,
00547                                  solver_params.first);
00548 
00549         totalrval.first  += rval.first;
00550         totalrval.second += rval.second;
00551       }
00552 
00553   this->release_linear_solver(linear_solver);
00554 
00555   for (unsigned int i=0; i != this->qoi.size(); ++i)
00556     if (qoi_indices.has_index(i))
00557       delete temprhs[i];
00558 
00559   // The linear solver may not have fit our constraints exactly
00560 #ifdef LIBMESH_ENABLE_CONSTRAINTS
00561   for (unsigned int i=0; i != this->qoi.size(); ++i)
00562     if (qoi_indices.has_index(i))
00563       this->get_dof_map().enforce_constraints_exactly
00564         (*this, &this->get_weighted_sensitivity_adjoint_solution(i),
00565          /* homogeneous = */ true);
00566 #endif
00567 
00568   // Stop logging the nonlinear solve
00569   STOP_LOG("weighted_sensitivity_adjoint_solve()", "ImplicitSystem");
00570 
00571   return totalrval;
00572 }

std::pair< unsigned int, Real > libMesh::ImplicitSystem::weighted_sensitivity_solve	(	const ParameterVector &	parameters,
		const ParameterVector &	weights
	)			`[virtual, inherited]`

Assembles & solves the linear system(s) (dR/du)*u_w = sum(w_p*-dR/dp), for those parameters p contained within parameters weighted by the values w_p found within weights.

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

Reimplemented from libMesh::System.

Definition at line 577 of file implicit_system.C.

References libMesh::System::add_weighted_sensitivity_solution(), libMesh::ImplicitSystem::assembly(), libMesh::NumericVector< T >::clone(), libMesh::SparseMatrix< T >::close(), libMesh::NumericVector< T >::close(), libMesh::ParameterVector::deep_copy(), libMesh::DofMap::enforce_constraints_exactly(), libMesh::System::get_dof_map(), libMesh::ImplicitSystem::get_linear_solve_parameters(), libMesh::ImplicitSystem::get_linear_solver(), libMesh::System::get_weighted_sensitivity_solution(), libMesh::ImplicitSystem::matrix, libMesh::Real, libMesh::ImplicitSystem::release_linear_solver(), libMesh::ExplicitSystem::rhs, libMesh::LinearSolver< T >::solve(), libMesh::TOLERANCE, and libMesh::ParameterVector::value_copy().

00579 {
00580   // Log how long the linear solve takes.
00581   START_LOG("weighted_sensitivity_solve()", "ImplicitSystem");
00582 
00583   // We currently get partial derivatives via central differencing
00584   const Real delta_p = TOLERANCE;
00585 
00586   // The forward system should now already be solved.
00587 
00588   // Now we're assembling a weighted sum of sensitivity systems:
00589   //
00590   // dR/du (u, v)(sum(w_l*u'_l)) = -sum_l(w_l*R'_l (u, v)) forall v
00591 
00592   // We'll assemble the rhs first, because the R' term will require
00593   // perturbing the system, and some applications may not be able to
00594   // assemble a perturbed residual without simultaneously constructing
00595   // a perturbed jacobian.
00596 
00597   // We approximate the _l partial derivatives via a central
00598   // differencing perturbation in the w_l direction:
00599   //
00600   // sum_l(w_l*v_l) ~= (v(p + dp*w_l*e_l) - v(p - dp*w_l*e_l))/(2*dp)
00601 
00602   ParameterVector oldparameters, parameterperturbation;
00603   parameters.deep_copy(oldparameters);
00604   weights.deep_copy(parameterperturbation);
00605   parameterperturbation *= delta_p;
00606   parameters += parameterperturbation;
00607 
00608   this->assembly(true, false);
00609   this->rhs->close();
00610 
00611   AutoPtr<NumericVector<Number> > temprhs = this->rhs->clone();
00612 
00613   oldparameters.value_copy(parameters);
00614   parameterperturbation *= -1.0;
00615   parameters += parameterperturbation;
00616 
00617   this->assembly(true, false);
00618   this->rhs->close();
00619 
00620   *temprhs -= *(this->rhs);
00621   *temprhs /= (2.0*delta_p);
00622 
00623   // Finally, assemble the jacobian at the non-perturbed parameter
00624   // values
00625   oldparameters.value_copy(parameters);
00626 
00627   // Build the Jacobian
00628   this->assembly(false, true);
00629   this->matrix->close();
00630 
00631   // The weighted sensitivity problem is linear
00632   LinearSolver<Number> *linear_solver = this->get_linear_solver();
00633 
00634   std::pair<unsigned int, Real> solver_params =
00635     this->get_linear_solve_parameters();
00636 
00637   const std::pair<unsigned int, Real> rval =
00638     linear_solver->solve (*matrix, this->add_weighted_sensitivity_solution(),
00639                           *temprhs,
00640                           solver_params.second,
00641                           solver_params.first);
00642 
00643   this->release_linear_solver(linear_solver);
00644 
00645   // The linear solver may not have fit our constraints exactly
00646 #ifdef LIBMESH_ENABLE_CONSTRAINTS
00647   this->get_dof_map().enforce_constraints_exactly
00648     (*this, &this->get_weighted_sensitivity_solution(),
00649      /* homogeneous = */ true);
00650 #endif
00651 
00652   // Stop logging the nonlinear solve
00653   STOP_LOG("weighted_sensitivity_solve()", "ImplicitSystem");
00654 
00655   return rval;
00656 }

void libMesh::System::write_header	(	Xdr &	io,
		const std::string &	version,
		const bool	write_additional_data
	)			const `[inherited]`

Writes the basic data header for this System.

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

for this system

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

for each variable in the system

6.) The name of the variable (string)

6.1.) subdomain where the variable lives

7.) Combined in an FEType:

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

end variable loop

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

for each additional vector in the system object

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

end system

Definition at line 1240 of file system_io.C.

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

01243 {
01277   libmesh_assert (io.writing());
01278 
01279 
01280   // Only write the header information
01281   // if we are processor 0.
01282   if (this->get_mesh().processor_id() != 0)
01283     return;
01284 
01285   std::string comment;
01286   char buf[80];
01287 
01288   // 5.)
01289   // Write the number of variables in the system
01290 
01291   {
01292     // set up the comment
01293     comment = "# No. of Variables in System \"";
01294     comment += this->name();
01295     comment += "\"";
01296 
01297     unsigned int nv = this->n_vars();
01298     io.data (nv, comment.c_str());
01299   }
01300 
01301 
01302   for (unsigned int var=0; var<this->n_vars(); var++)
01303     {
01304       // 6.)
01305       // Write the name of the var-th variable
01306       {
01307         // set up the comment
01308         comment  = "#   Name, Variable No. ";
01309         std::sprintf(buf, "%d", var);
01310         comment += buf;
01311         comment += ", System \"";
01312         comment += this->name();
01313         comment += "\"";
01314 
01315         std::string var_name = this->variable_name(var);
01316         io.data (var_name, comment.c_str());
01317       }
01318 
01319       // 6.1.) Variable subdomains
01320       {
01321         // set up the comment
01322                 comment  = "#     Subdomains, Variable \"";
01323                 std::sprintf(buf, "%s", this->variable_name(var).c_str());
01324                 comment += buf;
01325                 comment += "\", System \"";
01326                 comment += this->name();
01327                 comment += "\"";
01328 
01329                 const std::set<subdomain_id_type> & domains = this->variable(var).active_subdomains();
01330                 std::vector<subdomain_id_type> domain_array;
01331                 domain_array.assign(domains.begin(), domains.end());
01332                 io.data (domain_array, comment.c_str());
01333       }
01334 
01335       // 7.)
01336       // Write the approximation order of the var-th variable
01337       // in this system
01338       {
01339         // set up the comment
01340         comment = "#     Approximation Order, Variable \"";
01341         std::sprintf(buf, "%s", this->variable_name(var).c_str());
01342         comment += buf;
01343         comment += "\", System \"";
01344         comment += this->name();
01345         comment += "\"";
01346 
01347         int order = static_cast<int>(this->variable_type(var).order);
01348         io.data (order, comment.c_str());
01349       }
01350 
01351 
01352 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
01353 
01354       // do the same for radial_order
01355       {
01356         comment = "#     Radial Approximation Order, Variable \"";
01357         std::sprintf(buf, "%s", this->variable_name(var).c_str());
01358         comment += buf;
01359         comment += "\", System \"";
01360         comment += this->name();
01361         comment += "\"";
01362 
01363         int rad_order = static_cast<int>(this->variable_type(var).radial_order);
01364         io.data (rad_order, comment.c_str());
01365       }
01366 
01367 #endif
01368 
01369       // Write the Finite Element type of the var-th variable
01370       // in this System
01371       {
01372         // set up the comment
01373         comment = "#     FE Family, Variable \"";
01374         std::sprintf(buf, "%s", this->variable_name(var).c_str());
01375         comment += buf;
01376         comment += "\", System \"";
01377         comment += this->name();
01378         comment += "\"";
01379 
01380         const FEType& type = this->variable_type(var);
01381         int fam = static_cast<int>(type.family);
01382         io.data (fam, comment.c_str());
01383 
01384 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
01385 
01386         comment = "#     Radial FE Family, Variable \"";
01387         std::sprintf(buf, "%s", this->variable_name(var).c_str());
01388         comment += buf;
01389         comment += "\", System \"";
01390         comment += this->name();
01391         comment += "\"";
01392 
01393         int radial_fam = static_cast<int>(type.radial_family);
01394         io.data (radial_fam, comment.c_str());
01395 
01396         comment = "#     Infinite Mapping Type, Variable \"";
01397         std::sprintf(buf, "%s", this->variable_name(var).c_str());
01398         comment += buf;
01399         comment += "\", System \"";
01400         comment += this->name();
01401         comment += "\"";
01402 
01403         int i_map = static_cast<int>(type.inf_map);
01404         io.data (i_map, comment.c_str());
01405 #endif
01406       }
01407     } // end of the variable loop
01408 
01409   // 8.)
01410   // Write the number of additional vectors in the System.
01411   // If write_additional_data==false, then write zero for
01412   // the number of additional vectors.
01413   {
01414     {
01415       // set up the comment
01416       comment = "# No. of Additional Vectors, System \"";
01417       comment += this->name();
01418       comment += "\"";
01419 
01420       unsigned int nvecs = write_additional_data ? this->n_vectors () : 0;
01421       io.data (nvecs, comment.c_str());
01422     }
01423 
01424     if (write_additional_data)
01425       {
01426         std::map<std::string, NumericVector<Number>* >::const_iterator
01427           vec_pos = this->_vectors.begin();
01428         unsigned int cnt=0;
01429 
01430         for (; vec_pos != this->_vectors.end(); ++vec_pos)
01431           {
01432             // 9.)
01433             // write the name of the cnt-th additional vector
01434             comment =  "# Name of ";
01435             std::sprintf(buf, "%d", cnt++);
01436             comment += buf;
01437             comment += "th vector";
01438             std::string vec_name = vec_pos->first;
01439 
01440             io.data (vec_name, comment.c_str());
01441           }
01442       }
01443   }
01444 }

void libMesh::System::write_parallel_data	(	Xdr &	io,
		const bool	write_additional_data
	)			const `[inherited]`

Writes additional data, namely vectors, for this System. This method may safely be called on a distributed-memory mesh. This method will create an individual file for each processor in the simulation where the local solution components for that processor will be stored.

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

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

for each additional vector in the object

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

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

Definition at line 1448 of file system_io.C.

References libMesh::System::_vectors, libMesh::Xdr::data(), libMesh::FEType::family, libMesh::System::get_dof_map(), libMesh::System::get_mesh(), libMesh::DofObject::invalid_id, libMesh::n_processors(), libMesh::System::n_vars(), libMesh::System::name(), libMesh::System::number(), libMesh::processor_id(), libMeshEnums::SCALAR, libMesh::DofMap::SCALAR_dof_indices(), libMesh::System::solution, libMesh::Variable::type(), libMesh::System::variable(), and libMesh::Xdr::writing().

01450 {
01470   // PerfLog pl("IO Performance",false);
01471   // pl.push("write_parallel_data");
01472   std::size_t total_written_size = 0;
01473   
01474   std::string comment;
01475 
01476   libmesh_assert (io.writing());
01477 
01478   std::vector<Number> io_buffer; io_buffer.reserve(this->solution->local_size());
01479 
01480   // build the ordered nodes and element maps.
01481   // when writing/reading parallel files we need to iterate
01482   // over our nodes/elements in order of increasing global id().
01483   // however, this is not guaranteed to be ordering we obtain
01484   // by using the node_iterators/element_iterators directly.
01485   // so build a set, sorted by id(), that provides the ordering.
01486   // further, for memory economy build the set but then transfer
01487   // its contents to vectors, which will be sorted.
01488   std::vector<const DofObject*> ordered_nodes, ordered_elements;
01489   {
01490     std::set<const DofObject*, CompareDofObjectsByID>
01491       ordered_nodes_set (this->get_mesh().local_nodes_begin(),
01492                          this->get_mesh().local_nodes_end());
01493 
01494       ordered_nodes.insert(ordered_nodes.end(),
01495                            ordered_nodes_set.begin(),
01496                            ordered_nodes_set.end());
01497   }
01498   {
01499     std::set<const DofObject*, CompareDofObjectsByID>
01500       ordered_elements_set (this->get_mesh().local_elements_begin(),
01501                             this->get_mesh().local_elements_end());
01502 
01503       ordered_elements.insert(ordered_elements.end(),
01504                               ordered_elements_set.begin(),
01505                               ordered_elements_set.end());
01506   }
01507 
01508   const unsigned int sys_num = this->number();
01509   const unsigned int nv      = this->n_vars();
01510 
01511   // Loop over each non-SCALAR variable and each node, and write out the value.
01512   for (unsigned int var=0; var<nv; var++)
01513     if (this->variable(var).type().family != SCALAR)
01514     {
01515       // First write the node DOF values
01516       for (std::vector<const DofObject*>::const_iterator
01517              it = ordered_nodes.begin(); it != ordered_nodes.end(); ++it)
01518         for (unsigned int comp=0; comp<(*it)->n_comp(sys_num, var); comp++)
01519           {
01520             //libMesh::out << "(*it)->id()=" << (*it)->id() << std::endl;
01521             libmesh_assert_not_equal_to ((*it)->dof_number(sys_num, var, comp),
01522                                          DofObject::invalid_id);
01523 
01524             io_buffer.push_back((*this->solution)((*it)->dof_number(sys_num, var, comp)));
01525           }
01526 
01527       // Then write the element DOF values
01528       for (std::vector<const DofObject*>::const_iterator
01529              it = ordered_elements.begin(); it != ordered_elements.end(); ++it)
01530         for (unsigned int comp=0; comp<(*it)->n_comp(sys_num, var); comp++)
01531           {
01532             libmesh_assert_not_equal_to ((*it)->dof_number(sys_num, var, comp),
01533                                          DofObject::invalid_id);
01534 
01535             io_buffer.push_back((*this->solution)((*it)->dof_number(sys_num, var, comp)));
01536           }
01537     }
01538 
01539   // Finally, write the SCALAR data on the last processor
01540   for (unsigned int var=0; var<this->n_vars(); var++)
01541     if(this->variable(var).type().family == SCALAR)
01542       {
01543         if (libMesh::processor_id() == (libMesh::n_processors()-1))
01544           {
01545             const DofMap& dof_map = this->get_dof_map();
01546             std::vector<dof_id_type> SCALAR_dofs;
01547             dof_map.SCALAR_dof_indices(SCALAR_dofs, var);
01548 
01549             for(unsigned int i=0; i<SCALAR_dofs.size(); i++)
01550               {
01551                 io_buffer.push_back( (*this->solution)(SCALAR_dofs[i]) );
01552               }
01553           }
01554       }
01555 
01556   // 9.)
01557   //
01558   // Actually write the reordered solution vector
01559   // for the ith system to disk
01560 
01561   // set up the comment
01562   {
01563     comment = "# System \"";
01564     comment += this->name();
01565     comment += "\" Solution Vector";
01566   }
01567 
01568   io.data (io_buffer, comment.c_str());
01569 
01570   total_written_size += io_buffer.size();
01571 
01572   // Only write additional vectors if wanted
01573   if (write_additional_data)
01574     {
01575       std::map<std::string, NumericVector<Number>* >::const_iterator
01576         pos = _vectors.begin();
01577 
01578       for(; pos != this->_vectors.end(); ++pos)
01579         {
01580           io_buffer.clear(); io_buffer.reserve( pos->second->local_size());
01581 
01582           // Loop over each non-SCALAR variable and each node, and write out the value.
01583           for (unsigned int var=0; var<nv; var++)
01584             if(this->variable(var).type().family != SCALAR)
01585             {
01586               // First write the node DOF values
01587               for (std::vector<const DofObject*>::const_iterator
01588                      it = ordered_nodes.begin(); it != ordered_nodes.end(); ++it)
01589                 for (unsigned int comp=0; comp<(*it)->n_comp(sys_num, var); comp++)
01590                   {
01591                     libmesh_assert_not_equal_to ((*it)->dof_number(sys_num, var, comp),
01592                                                  DofObject::invalid_id);
01593 
01594                     io_buffer.push_back((*pos->second)((*it)->dof_number(sys_num, var, comp)));
01595                   }
01596 
01597               // Then write the element DOF values
01598               for (std::vector<const DofObject*>::const_iterator
01599                      it = ordered_elements.begin(); it != ordered_elements.end(); ++it)
01600                 for (unsigned int comp=0; comp<(*it)->n_comp(sys_num, var); comp++)
01601                   {
01602                     libmesh_assert_not_equal_to ((*it)->dof_number(sys_num, var, comp),
01603                                                  DofObject::invalid_id);
01604 
01605                     io_buffer.push_back((*pos->second)((*it)->dof_number(sys_num, var, comp)));
01606                   }
01607             }
01608 
01609             // Finally, write the SCALAR data on the last processor
01610             for (unsigned int var=0; var<this->n_vars(); var++)
01611               if(this->variable(var).type().family == SCALAR)
01612                 {
01613                   if (libMesh::processor_id() == (libMesh::n_processors()-1))
01614                     {
01615                       const DofMap& dof_map = this->get_dof_map();
01616                       std::vector<dof_id_type> SCALAR_dofs;
01617                       dof_map.SCALAR_dof_indices(SCALAR_dofs, var);
01618 
01619                       for(unsigned int i=0; i<SCALAR_dofs.size(); i++)
01620                         {
01621                           io_buffer.push_back( (*pos->second)(SCALAR_dofs[i]) );
01622                         }
01623                     }
01624                 }
01625 
01626           // 10.)
01627           //
01628           // Actually write the reordered additional vector
01629           // for this system to disk
01630 
01631           // set up the comment
01632           {
01633             comment = "# System \"";
01634             comment += this->name();
01635             comment += "\" Additional Vector \"";
01636             comment += pos->first;
01637             comment += "\"";
01638           }
01639 
01640           io.data (io_buffer, comment.c_str());
01641 
01642           total_written_size += io_buffer.size();
01643         }
01644     }
01645 
01646   // const Real 
01647   //   dt   = pl.get_elapsed_time(),
01648   //   rate = total_written_size*sizeof(Number)/dt; 
01649 
01650   // std::cerr << "Write " << total_written_size << " \"Number\" values\n"
01651   //        << " Elapsed time = " << dt << '\n'
01652   //        << " Rate = " << rate/1.e6 << "(MB/sec)\n\n";
01653   
01654   // pl.pop("write_parallel_data");
01655 }

void libMesh::System::write_serialized_data	(	Xdr &	io,
		const bool	write_additional_data = `true`
	)			const `[inherited]`

Writes additional data, namely vectors, for this System. This method may safely be called on a distributed-memory mesh.

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

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

for each additional vector in the object

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

Definition at line 1659 of file system_io.C.

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

01661 {
01675   parallel_only();
01676   std::string comment;
01677 
01678   // PerfLog pl("IO Performance",false);
01679   // pl.push("write_serialized_data");
01680   std::size_t total_written_size = 0;
01681   
01682   total_written_size +=
01683     this->write_serialized_vector(io, *this->solution);
01684 
01685   // set up the comment
01686   if (libMesh::processor_id() == 0)
01687     {
01688       comment = "# System \"";
01689       comment += this->name();
01690       comment += "\" Solution Vector";
01691 
01692       io.comment (comment);
01693     }
01694 
01695   // Only write additional vectors if wanted
01696   if (write_additional_data)
01697     {
01698       std::map<std::string, NumericVector<Number>* >::const_iterator
01699         pos = _vectors.begin();
01700 
01701       for(; pos != this->_vectors.end(); ++pos)
01702         {
01703           total_written_size +=
01704             this->write_serialized_vector(io, *pos->second);
01705 
01706           // set up the comment
01707           if (libMesh::processor_id() == 0)
01708             {
01709               comment = "# System \"";
01710               comment += this->name();
01711               comment += "\" Additional Vector \"";
01712               comment += pos->first;
01713               comment += "\"";
01714               io.comment (comment);
01715             }
01716         }
01717     }
01718 
01719   // const Real 
01720   //   dt   = pl.get_elapsed_time(),
01721   //   rate = total_written_size*sizeof(Number)/dt; 
01722 
01723   // std::cout << "Write " << total_written_size << " \"Number\" values\n"
01724   //        << " Elapsed time = " << dt << '\n'
01725   //        << " Rate = " << rate/1.e6 << "(MB/sec)\n\n";
01726   
01727   // pl.pop("write_serialized_data");
01728 
01729   
01730 
01731   
01732   // // test the new method
01733   // {
01734   //   std::vector<std::string> names;
01735   //   std::vector<NumericVector<Number>*> vectors_to_write;
01736 
01737   //   names.push_back("Solution Vector");
01738   //   vectors_to_write.push_back(this->solution.get());
01739     
01740   //   // Only write additional vectors if wanted
01741   //   if (write_additional_data)
01742   //     {
01743   //    std::map<std::string, NumericVector<Number>* >::const_iterator
01744   //      pos = _vectors.begin();
01745         
01746   //    for(; pos != this->_vectors.end(); ++pos)
01747   //      {
01748   //        names.push_back("Additional Vector " + pos->first);
01749   //        vectors_to_write.push_back(pos->second);
01750   //      }
01751   //     }
01752 
01753   //   total_written_size =
01754   //     this->write_serialized_vectors (io, names, vectors_to_write);
01755 
01756   //   const Real 
01757   //     dt2   = pl.get_elapsed_time(),
01758   //     rate2 = total_written_size*sizeof(Number)/(dt2-dt);
01759     
01760   //   std::cout << "Write (new) " << total_written_size << " \"Number\" values\n"
01761   //          << " Elapsed time = " << (dt2-dt) << '\n'
01762   //          << " Rate = " << rate2/1.e6 << "(MB/sec)\n\n";
01763 
01764   // }
01765 }

dof_id_type libMesh::System::write_serialized_vectors	(	Xdr &	io,
		const std::vector< const NumericVector< Number > * > &	vectors
	)			const `[inherited]`

Serialize & write a number of identically distributed vectors. This method allows for optimization for the multiple vector case by only communicating the metadata once.

Definition at line 2256 of file system_io.C.

References libMesh::Xdr::data(), libMesh::FEType::family, libMesh::System::get_mesh(), libMesh::MeshBase::n_elem(), libMesh::MeshTools::n_elem(), libMesh::MeshBase::n_nodes(), n_nodes, libMesh::System::n_vars(), libMesh::processor_id(), libMeshEnums::SCALAR, libMesh::Variable::type(), libMesh::System::variable(), libMesh::System::write_SCALAR_dofs(), libMesh::System::write_serialized_blocked_dof_objects(), and libMesh::Xdr::writing().

02258 {
02259   parallel_only();
02260 
02261   libmesh_assert (io.writing());
02262 
02263   // Cache these - they are not free!
02264   const dof_id_type
02265     n_nodes       = this->get_mesh().n_nodes(),
02266     n_elem        = this->get_mesh().n_elem();  
02267 
02268   dof_id_type written_length = 0.;
02269   
02270   if (libMesh::processor_id() == 0) 
02271     {
02272       unsigned int
02273         n_vec    = libmesh_cast_int<unsigned int>(vectors.size());
02274       dof_id_type
02275         vec_size = vectors.empty() ? 0 : vectors[0]->size();
02276       // Set the number of vectors
02277       io.data(n_vec, "# number of vectors");
02278       // Set the buffer size
02279       io.data(vec_size, "# vector length");
02280     }
02281 
02282   //---------------------------------
02283   // Collect the values for all nodes
02284   written_length +=
02285     this->write_serialized_blocked_dof_objects (vectors,
02286                                                 n_nodes,
02287                                                 this->get_mesh().local_nodes_begin(),
02288                                                 this->get_mesh().local_nodes_end(),
02289                                                 io);
02290 
02291   //------------------------------------
02292   // Collect the values for all elements
02293   written_length +=
02294     this->write_serialized_blocked_dof_objects (vectors,
02295                                                 n_elem,
02296                                                 this->get_mesh().local_elements_begin(),
02297                                                 this->get_mesh().local_elements_end(),
02298                                                 io);
02299 
02300   //-------------------------------------------
02301   // Finally loop over all the SCALAR variables
02302   for (unsigned int vec=0; vec<vectors.size(); vec++)
02303     for (unsigned int var=0; var<this->n_vars(); var++)
02304       if(this->variable(var).type().family == SCALAR)
02305         {
02306           libmesh_assert_not_equal_to (vectors[vec], 0);
02307           
02308           written_length +=
02309             this->write_SCALAR_dofs (*vectors[vec], var, io);
02310         }
02311       
02312   return written_length;
02313 }

void libMesh::System::zero_variable	(	NumericVector< Number > &	v,
		unsigned int	var_num
	)			const `[inherited]`

Zeroes all dofs in v that correspond to variable number var_num.

Definition at line 1314 of file system.C.

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

01315 {
01316   /* Make sure the call makes sense.  */
01317   libmesh_assert_less (var_num, this->n_vars());
01318 
01319   /* Get a reference to the mesh.  */
01320   const MeshBase& mesh = this->get_mesh();
01321 
01322   /* Check which system we are.  */
01323   const unsigned int sys_num = this->number();
01324 
01325   /* Loop over nodes.  */
01326   {
01327     MeshBase::const_node_iterator it = mesh.local_nodes_begin();
01328     const MeshBase::const_node_iterator end_it = mesh.local_nodes_end();
01329     for ( ; it != end_it; ++it)
01330       {
01331         const Node* node = *it;
01332         unsigned int n_comp = node->n_comp(sys_num,var_num);
01333         for(unsigned int i=0; i<n_comp; i++)
01334           {
01335             const dof_id_type index = node->dof_number(sys_num,var_num,i);
01336             v.set(index,0.0);
01337           }
01338       }
01339   }
01340 
01341   /* Loop over elements.  */
01342   {
01343     MeshBase::const_element_iterator it = mesh.active_local_elements_begin();
01344     const MeshBase::const_element_iterator end_it = mesh.active_local_elements_end();
01345     for ( ; it != end_it; ++it)
01346       {
01347         const Elem* elem = *it;
01348         unsigned int n_comp = elem->n_comp(sys_num,var_num);
01349         for(unsigned int i=0; i<n_comp; i++)
01350           {
01351             const dof_id_type index = elem->dof_number(sys_num,var_num,i);
01352             v.set(index,0.0);
01353           }
01354       }
01355   }
01356 }

Member Data Documentation

ReferenceCounter::Counts libMesh::ReferenceCounter::_counts [static, protected, inherited]

Actually holds the data.

Definition at line 118 of file reference_counter.h.

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

bool libMesh::ReferenceCounter::_enable_print_counter = true [static, protected, inherited]

Flag to control whether reference count information is printed when print_info is called.

Definition at line 137 of file reference_counter.h.

Referenced by libMesh::ReferenceCounter::disable_print_counter_info(), libMesh::ReferenceCounter::enable_print_counter_info(), and libMesh::ReferenceCounter::print_info().

Real libMesh::LinearImplicitSystem::_final_linear_residual [protected]

The final residual for the linear system Ax=b.

Definition at line 195 of file linear_implicit_system.h.

Referenced by final_linear_residual(), solve(), and libMesh::FrequencySystem::solve().

Threads::spin_mutex libMesh::ReferenceCounter::_mutex [static, protected, inherited]

Mutual exclusion object to enable thread-safe reference counting.

Definition at line 131 of file reference_counter.h.

unsigned int libMesh::LinearImplicitSystem::_n_linear_iterations [protected]

The number of linear iterations required to solve the linear system Ax=b.

Definition at line 190 of file linear_implicit_system.h.

Referenced by n_linear_iterations(), solve(), and libMesh::FrequencySystem::solve().

Threads::atomic< unsigned int > libMesh::ReferenceCounter::_n_objects [static, protected, inherited]

The number of objects. Print the reference count information when the number returns to 0.

Definition at line 126 of file reference_counter.h.

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

ShellMatrix<Number>* libMesh::LinearImplicitSystem::_shell_matrix [protected]

User supplies shell matrix or NULL if no shell matrix is used.

Definition at line 200 of file linear_implicit_system.h.

Referenced by attach_shell_matrix(), get_shell_matrix(), and solve().

const SystemSubset* libMesh::LinearImplicitSystem::_subset [protected]

The current subset on which to solve (or NULL if none).

Definition at line 205 of file linear_implicit_system.h.

Referenced by restrict_solve_to(), and solve().

SubsetSolveMode libMesh::LinearImplicitSystem::_subset_solve_mode [protected]

If restrict-solve-to-subset mode is active, this member decides what happens with the dofs outside the subset.

Definition at line 211 of file linear_implicit_system.h.

Referenced by restrict_solve_to(), and solve().

bool libMesh::System::assemble_before_solve [inherited]

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

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

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

Definition at line 1371 of file system.h.

Referenced by libMesh::ImplicitSystem::adjoint_solve(), libMesh::System::disable_cache(), libMesh::ImplicitSystem::disable_cache(), libMesh::ImplicitSystem::sensitivity_solve(), solve(), libMesh::EigenSystem::solve(), and libMesh::CondensedEigenSystem::solve().

AutoPtr<NumericVector<Number> > libMesh::System::current_local_solution [inherited]

All the values I need to compute my contribution to the simulation at hand. Think of this as the current solution with any ghost values needed from other processors. This vector is necessarily larger than the solution vector in the case of a parallel simulation. The update() member is used to synchronize the contents of the solution and current_local_solution vectors.

Definition at line 1430 of file system.h.

Referenced by libMesh::NonlinearImplicitSystem::assembly(), libMesh::System::clear(), libMesh::Problem_Interface::computeF(), libMesh::Problem_Interface::computeJacobian(), libMesh::Problem_Interface::computePreconditioner(), libMesh::System::current_solution(), DMFunction_libMesh(), DMJacobian_libMesh(), libMesh::ExactErrorEstimator::estimate_error(), libMesh::AdjointRefinementEstimator::estimate_error(), libMesh::System::init_data(), libMesh::System::project_solution(), libMesh::System::re_update(), libMesh::System::reinit(), libMesh::System::restrict_vectors(), and libMesh::System::update().

int libMesh::System::extra_quadrature_order [inherited]

A member int that can be employed to indicate increased or reduced quadrature order.

Note for FEMSystem users: By default, when calling the user-defined residual functions, the FEMSystem will first set up an appropriate FEType::default_quadrature_rule() object for performing the integration. This rule will integrate elements of order up to 2*p+1 exactly (where p is the sum of the base FEType and local p refinement levels), but if additional (or reduced) quadrature accuracy is desired then this extra_quadrature_order (default 0) will be added.

Definition at line 1403 of file system.h.

Referenced by libMesh::FEMContext::FEMContext().

AutoPtr<LinearSolver<Number> > libMesh::LinearImplicitSystem::linear_solver

The LinearSolver defines the interface used to solve the linear_implicit system. This class handles all the details of interfacing with various linear algebra packages like PETSc or LASPACK.

Definition at line 149 of file linear_implicit_system.h.

Referenced by clear(), get_linear_solver(), init_data(), reinit(), solve(), and libMesh::FrequencySystem::solve().

SparseMatrix<Number>* libMesh::ImplicitSystem::matrix [inherited]

The system matrix. Implicit systems are characterized by the need to solve the linear system Ax=b. This is the system matrix A.

Definition at line 335 of file implicit_system.h.

Referenced by libMesh::__libmesh_petsc_diff_solver_jacobian(), libMesh::ImplicitSystem::add_system_matrix(), libMesh::ImplicitSystem::adjoint_solve(), libMesh::ImplicitSystem::assemble(), libMesh::NonlinearImplicitSystem::assembly(), assembly(), libMesh::FEMSystem::assembly(), libMesh::ImplicitSystem::clear(), libMesh::NewmarkSystem::compute_matrix(), libMesh::ContinuationSystem::continuation_solve(), DMCreateMatrix_libMesh(), DMJacobian_libMesh(), libMesh::ImplicitSystem::init_matrices(), libMesh::ImplicitSystem::sensitivity_solve(), libMesh::NoxNonlinearSolver< T >::solve(), libMesh::PetscDiffSolver::solve(), libMesh::NonlinearImplicitSystem::solve(), libMesh::NewtonSolver::solve(), solve(), libMesh::FrequencySystem::solve(), libMesh::EigenTimeSolver::solve(), libMesh::ContinuationSystem::solve_tangent(), libMesh::ImplicitSystem::weighted_sensitivity_adjoint_solve(), and libMesh::ImplicitSystem::weighted_sensitivity_solve().

std::vector<Number> libMesh::System::qoi [inherited]

Values of the quantities of interest. This vector needs to be both resized and filled by the user before any quantity of interest assembly is done and before any sensitivities are calculated.

Definition at line 1449 of file system.h.

Referenced by libMesh::ImplicitSystem::adjoint_solve(), libMesh::SensitivityData::allocate_data(), libMesh::SensitivityData::allocate_hessian_data(), libMesh::FEMSystem::assemble_qoi(), libMesh::ExplicitSystem::assemble_qoi(), libMesh::FEMSystem::assemble_qoi_derivative(), libMesh::ExplicitSystem::assemble_qoi_derivative(), libMesh::DifferentiableSystem::attach_qoi(), libMesh::DiffContext::DiffContext(), libMesh::AdjointResidualErrorEstimator::estimate_error(), libMesh::AdjointRefinementEstimator::estimate_error(), libMesh::FEMContext::pre_fe_reinit(), libMesh::QoISet::size(), and libMesh::ImplicitSystem::weighted_sensitivity_adjoint_solve().