FEMSystem Class Reference

#include <fem_system.h>

Inheritance diagram for FEMSystem:

List of all members.

Public Types

typedef FEMSystem sys_type

typedef DifferentiableSystem 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

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

virtual ~FEMSystem ()

virtual void clear ()

virtual void assembly (bool get_residual, bool get_jacobian)

virtual void solve ()

virtual void time_evolving (unsigned int var)

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

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

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

void mesh_position_get ()

void mesh_position_set ()

virtual bool eulerian_residual (bool request_jacobian, DiffContext &context)

virtual bool mass_residual (bool request_jacobian, DiffContext &context)

virtual AutoPtr< DiffContext > build_context ()

virtual void init_context (DiffContext &)

virtual void postprocess ()

virtual void assemble_qoi (const QoISet &indices=QoISet())

virtual void assemble_qoi_derivative (const QoISet &indices=QoISet())

virtual void reinit ()

virtual void assemble ()

virtual LinearSolver< Number > * get_linear_solver () const

virtual std::pair< unsigned
int, Real > get_linear_solve_parameters () const

virtual void release_linear_solver (LinearSolver< Number > *) const

virtual bool element_time_derivative (bool request_jacobian, DiffContext &)

virtual bool element_constraint (bool request_jacobian, DiffContext &)

virtual bool side_time_derivative (bool request_jacobian, DiffContext &)

virtual bool side_constraint (bool request_jacobian, DiffContext &)

virtual void element_postprocess (DiffContext &)

virtual void side_postprocess (DiffContext &)

virtual void element_qoi (DiffContext &, const QoISet &)

virtual void element_qoi_derivative (DiffContext &, const QoISet &)

virtual void side_qoi (DiffContext &, const QoISet &)

virtual void side_qoi_derivative (DiffContext &, const QoISet &)

virtual bool side_mass_residual (bool request_jacobian, DiffContext &)

sys_type & system ()

virtual std::string system_type () 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)

unsigned int n_matrices () const

void init ()

virtual void update ()

virtual void qoi_parameter_sensitivity (const QoISet &qoi_indices, const ParameterVector &parameters, SensitivityData &sensitivities)

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

const std::string & name () const

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

void project_vector (Number fptr(const Point &p, const Parameters &parameters, const std::string &sys_name, const std::string &unknown_name), Gradient gptr(const Point &p, const Parameters &parameters, const std::string &sys_name, const std::string &unknown_name), Parameters &parameters, NumericVector< Number > &new_vector) const

unsigned int number () const

void update_global_solution (std::vector< Number > &global_soln) const

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

const MeshBase & get_mesh () const

MeshBase & get_mesh ()

const DofMap & get_dof_map () const

DofMap & get_dof_map ()

const EquationSystems & get_equation_systems () const

EquationSystems & get_equation_systems ()

bool active () const

void activate ()

void deactivate ()

vectors_iterator vectors_begin ()

const_vectors_iterator vectors_begin () const

vectors_iterator vectors_end ()

const_vectors_iterator vectors_end () const

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

bool & project_solution_on_reinit (void)

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

const NumericVector< Number > * request_vector (const std::string &vec_name) const

NumericVector< Number > * request_vector (const std::string &vec_name)

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

NumericVector< Number > & add_weighted_sensitivity_solution ()

NumericVector< Number > & get_weighted_sensitivity_solution ()

const NumericVector< Number > & get_weighted_sensitivity_solution () const

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

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

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

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

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

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

unsigned int n_vectors () const

unsigned int n_vars () const

unsigned int n_dofs () const

unsigned int n_active_dofs () const

unsigned int n_constrained_dofs () const

unsigned int n_local_dofs () const

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

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

const Variable & variable (unsigned int var) const

bool has_variable (const std::string &var) const

const std::string & variable_name (const unsigned int i) const

unsigned short int variable_number (const std::string &var) const

const FEType & variable_type (const unsigned int i) const

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

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

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

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

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

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

void read_parallel_data (Xdr &io, const bool read_additional_data)

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

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

void write_parallel_data (Xdr &io, const bool write_additional_data) const

std::string get_info () const

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

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

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

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

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

virtual void user_initialization ()

virtual void user_assembly ()

virtual void user_constrain ()

virtual void user_QOI (const QoISet &qoi_indices)

virtual void user_QOI_derivative (const QoISet &qoi_indices)

virtual void re_update ()

virtual void restrict_vectors ()

virtual void prolong_vectors ()

Number current_solution (const unsigned int global_dof_number) const

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

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

Static Public Member Functions

static std::string get_info ()

static void print_info ()

static unsigned int n_objects ()

Public Attributes

bool fe_reinit_during_postprocess

int extra_quadrature_order

Real numerical_jacobian_h

Real verify_analytic_jacobians

bool compute_internal_sides

bool postprocess_sides

bool assemble_qoi_sides

AutoPtr< TimeSolver > time_solver

Real time

Real deltat

bool print_solution_norms

bool print_solutions

bool print_residual_norms

bool print_residuals

bool print_jacobian_norms

bool print_jacobians

bool print_element_jacobians

bool use_fixed_solution

SparseMatrix< Number > * matrix

NumericVector< Number > * rhs

bool assemble_before_solve

AutoPtr< NumericVector< Number > > solution

AutoPtr< NumericVector< Number > > current_local_solution

std::vector< Number > qoi

Protected Types

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

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

Protected Member Functions

virtual void init_data ()

void numerical_jacobian (TimeSolverResPtr res, FEMContext &context)

void numerical_elem_jacobian (FEMContext &context)

void numerical_side_jacobian (FEMContext &context)

virtual void init_matrices ()

void project_vector (NumericVector< Number > &) const

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

void increment_constructor_count (const std::string &name)

void increment_destructor_count (const std::string &name)

Protected Attributes

unsigned int _mesh_sys

unsigned int _mesh_x_var

unsigned int _mesh_y_var

unsigned int _mesh_z_var

std::vector< bool > _time_evolving

Static Protected Attributes

static Counts _counts

static Threads::atomic
< unsigned int > _n_objects

static Threads::spin_mutex _mutex

Detailed Description

This class provides a specific system class. It aims at nonlinear implicit systems, requiring only a cell residual calculation from the user. Note that still additional vectors/matrices may be added, as offered in the class ExplicitSystem.

This class is part of the new DifferentiableSystem framework, which is still experimental. Users of this framework should beware of bugs and future API changes.

Author:: Roy H. Stogner 2006

Definition at line 52 of file fem_system.h.

Member Typedef Documentation

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

Definition at line 268 of file implicit_system.h.

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

Definition at line 411 of file system.h.

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

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

Definition at line 105 of file reference_counter.h.

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

Matrix iterator typedefs.

Definition at line 267 of file implicit_system.h.

typedef DifferentiableSystem FEMSystem::Parent

The type of the parent.

Reimplemented from DifferentiableSystem.

Reimplemented in ContinuationSystem.

Definition at line 77 of file fem_system.h.

typedef FEMSystem FEMSystem::sys_type

The type of system.

Reimplemented from DifferentiableSystem.

Reimplemented in ContinuationSystem.

Definition at line 72 of file fem_system.h.

typedef bool(TimeSolver::* FEMSystem::TimeSolverResPtr)(bool, DiffContext &) [protected]

Syntax sugar to make numerical_jacobian() declaration easier.

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

Vector iterator typedefs.

Definition at line 410 of file system.h.

Constructor & Destructor Documentation

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

Constructor. Optionally initializes required data structures.

Definition at line 41 of file fem_system.C.

00044   : Parent(es, name, number),
00045     fe_reinit_during_postprocess(true),
00046     extra_quadrature_order(0),
00047     numerical_jacobian_h(1.e-6),
00048     verify_analytic_jacobians(0.0),
00049     _mesh_sys(libMesh::invalid_uint),
00050     _mesh_x_var(libMesh::invalid_uint),
00051     _mesh_y_var(libMesh::invalid_uint),
00052     _mesh_z_var(libMesh::invalid_uint)
00053 {
00054 }

FEMSystem::~FEMSystem ( ) [virtual]

Destructor.

Definition at line 57 of file fem_system.C.

References clear().

00058 {
00059   this->clear();
00060 }

Member Function Documentation

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

Activates the system. Only active systems are solved.

Definition at line 1380 of file system.h.

References System::_active.

01381 {
01382   _active = true;
01383 }

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

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

Definition at line 1372 of file system.h.

References System::_active.

01373 {
01374   return _active;  
01375 }

NumericVector< Number > & System::add_adjoint_rhs ( unsigned int i = 0 ) [inherited]

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

Definition at line 824 of file system.C.

References System::add_vector().

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

00825 {
00826   OStringStream adjoint_rhs_name;
00827   adjoint_rhs_name << "adjoint_rhs" << i;
00828 
00829   return this->add_vector(adjoint_rhs_name.str());
00830 }

NumericVector< Number > & System::add_adjoint_solution ( unsigned int i = 0 ) [inherited]

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

Definition at line 764 of file system.C.

References System::add_vector().

Referenced by ImplicitSystem::adjoint_solve().

00765 {
00766   OStringStream adjoint_name;
00767   adjoint_name << "adjoint_solution" << i;
00768 
00769   return this->add_vector(adjoint_name.str());
00770 }

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

Adds the additional matrix mat_name to this system. Only allowed prior to assemble(). All additional matrices have the same sparsity pattern as the matrix used during solution. When not System but the user wants to initialize the mayor matrix, then all the additional matrices, if existent, have to be initialized by the user, too.

Definition at line 175 of file implicit_system.C.

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

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

00176 {
00177   // only add matrices before initializing...
00178   if (!_can_add_matrices)
00179     {
00180       std::cerr << "ERROR: Too late.  Cannot add matrices to the system after initialization"
00181                 << std::endl
00182                 << " any more.  You should have done this earlier."
00183                 << std::endl;
00184       libmesh_error();
00185     }
00186 
00187   // Return the matrix if it is already there.
00188   if (this->have_matrix(mat_name))
00189     return *(_matrices[mat_name]);
00190 
00191   // Otherwise build the matrix and return it.
00192   SparseMatrix<Number>* buf = SparseMatrix<Number>::build().release();
00193   _matrices.insert (std::make_pair (mat_name, buf));
00194 
00195   return *buf;
00196 }

NumericVector< Number > & System::add_sensitivity_rhs ( unsigned int i = 0 ) [inherited]

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

Definition at line 854 of file system.C.

References System::add_vector().

Referenced by ImplicitSystem::assemble_residual_derivatives().

00855 {
00856   OStringStream sensitivity_rhs_name;
00857   sensitivity_rhs_name << "sensitivity_rhs" << i;
00858 
00859   return this->add_vector(sensitivity_rhs_name.str());
00860 }

NumericVector< Number > & System::add_sensitivity_solution ( unsigned int i = 0 ) [inherited]

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

Definition at line 714 of file system.C.

References System::add_vector().

00715 {
00716   OStringStream sensitivity_name;
00717   sensitivity_name << "sensitivity_solution" << i;
00718 
00719   return this->add_vector(sensitivity_name.str());
00720 }

unsigned int System::add_variable	(	const std::string &	var,
		const Order	order = `FIRST`,
		const FEFamily	family = `LAGRANGE`,
		const std::set< subdomain_id_type > *const	active_subdomains = `NULL`
	)			`[inherited]`

Adds the variable var to the list of variables for this system. Same as before, but assumes LAGRANGE as default value for FEType.family.

Definition at line 923 of file system.C.

References System::add_variable().

00927 {
00928   return this->add_variable(var, 
00929                             FEType(order, family), 
00930                             active_subdomains);
00931 }

unsigned int System::add_variable	(	const std::string &	var,
		const FEType &	type,
		const std::set< subdomain_id_type > *const	active_subdomains = `NULL`
	)			`[inherited]`

Adds the variable var to the list of variables for this system. Returns the index number for the new variable.

Definition at line 884 of file system.C.

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

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

00887 {  
00888   // Make sure the variable isn't there already
00889   // or if it is, that it's the type we want
00890   for (unsigned int v=0; v<this->n_vars(); v++)
00891     if (this->variable_name(v) == var)
00892       {
00893         if (this->variable_type(v) == type)
00894           return _variables[v].number();
00895 
00896         std::cerr << "ERROR: incompatible variable "
00897                   << var
00898                   << " has already been added for this system!"
00899                   << std::endl;
00900         libmesh_error();
00901       }
00902 
00903   const unsigned int curr_n_vars = this->n_vars();                                              
00904 
00905   // Add the variable to the list
00906   _variables.push_back((active_subdomains == NULL) ?
00907                        Variable(var, curr_n_vars, type) :
00908                        Variable(var, curr_n_vars, type, *active_subdomains));
00909 
00910   libmesh_assert ((curr_n_vars+1) == this->n_vars());
00911 
00912   _variable_numbers[var] = curr_n_vars;
00913 
00914   // Add the variable to the _dof_map
00915   _dof_map->add_variable (_variables.back());
00916 
00917   // Return the number of the new variable
00918   return curr_n_vars;
00919 }

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

Adds the additional vector vec_name to this system. Only allowed prior to init(). All the additional vectors are similarly distributed, like the solution, and inititialized to zero.

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

Definition at line 549 of file system.C.

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

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

00551 {
00552   // Return the vector if it is already there.
00553   if (this->have_vector(vec_name))
00554     return *(_vectors[vec_name]);
00555 
00556   // Otherwise build the vector
00557   NumericVector<Number>* buf = NumericVector<Number>::build().release();
00558   _vectors.insert (std::make_pair (vec_name, buf));
00559   _vector_projections.insert (std::make_pair (vec_name, projections));
00560 
00561   // Initialize it if necessary
00562   if (!_can_add_vectors)
00563     buf->init (this->n_dofs(), this->n_local_dofs(), false, PARALLEL);
00564 
00565   return *buf;
00566 }

NumericVector< Number > & System::add_weighted_sensitivity_adjoint_solution ( unsigned int i = 0 ) [inherited]

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

Definition at line 794 of file system.C.

References System::get_vector().

Referenced by ImplicitSystem::weighted_sensitivity_adjoint_solve().

00795 {
00796   OStringStream adjoint_name;
00797   adjoint_name << "weighted_sensitivity_adjoint_solution" << i;
00798 
00799   return this->get_vector(adjoint_name.str());
00800 }

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

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

Definition at line 743 of file system.C.

References System::add_vector().

Referenced by ImplicitSystem::weighted_sensitivity_solve().

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

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

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

Uses adjoint_solve() and the adjoint sensitivity method.

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

Reimplemented from System.

Definition at line 657 of file implicit_system.C.

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

00660 {
00661   const unsigned int Np = parameters.size();
00662   const unsigned int Nq = qoi.size();
00663 
00664   // We currently get partial derivatives via central differencing
00665   const Real delta_p = TOLERANCE;
00666 
00667   // An introduction to the problem:
00668   //
00669   // Residual R(u(p),p) = 0
00670   // partial R / partial u = J = system matrix
00671   //
00672   // This implies that:
00673   // d/dp(R) = 0
00674   // (partial R / partial p) + 
00675   // (partial R / partial u) * (partial u / partial p) = 0
00676 
00677   // We first do an adjoint solve:
00678   // J^T * z = (partial q / partial u)
00679 
00680   this->adjoint_solve(qoi_indices);
00681 
00682   // Get ready to fill in senstivities:
00683   sensitivities.allocate_data(qoi_indices, *this, parameters);
00684 
00685   // We use the identities:
00686   // dq/dp = (partial q / partial p) + (partial q / partial u) *
00687   //         (partial u / partial p)
00688   // dq/dp = (partial q / partial p) + (J^T * z) *
00689   //         (partial u / partial p)
00690   // dq/dp = (partial q / partial p) + z * J *
00691   //         (partial u / partial p)
00692  
00693   // Leading to our final formula:
00694   // dq/dp = (partial q / partial p) - z * (partial R / partial p)
00695 
00696   for (unsigned int j=0; j != Np; ++j)
00697     {
00698       // (partial q / partial p) ~= (q(p+dp)-q(p-dp))/(2*dp)
00699       // (partial R / partial p) ~= (rhs(p+dp) - rhs(p-dp))/(2*dp)
00700 
00701       Number old_parameter = *parameters[j];
00702       // Number old_qoi = this->qoi;
00703 
00704       *parameters[j] = old_parameter - delta_p;
00705       this->assemble_qoi(qoi_indices);
00706       std::vector<Number> qoi_minus = this->qoi;
00707 
00708       this->assembly(true, false);
00709       this->rhs->close();
00710       AutoPtr<NumericVector<Number> > partialR_partialp = this->rhs->clone();
00711       *partialR_partialp *= -1;
00712 
00713       *parameters[j] = old_parameter + delta_p;
00714       this->assemble_qoi(qoi_indices);
00715       std::vector<Number>& qoi_plus = this->qoi;
00716 
00717       std::vector<Number> partialq_partialp(Nq, 0);
00718       for (unsigned int i=0; i != Nq; ++i)
00719         if (qoi_indices.has_index(i))
00720           partialq_partialp[i] = (qoi_plus[i] - qoi_minus[i]) / (2.*delta_p);
00721 
00722       this->assembly(true, false);
00723       this->rhs->close();
00724       *partialR_partialp += *this->rhs;
00725       *partialR_partialp /= (2.*delta_p);
00726 
00727       // Don't leave the parameter changed
00728       *parameters[j] = old_parameter;
00729 
00730       for (unsigned int i=0; i != Nq; ++i)
00731         if (qoi_indices.has_index(i))
00732           sensitivities[i][j] = partialq_partialp[i] -
00733             partialR_partialp->dot(this->get_adjoint_solution(i));
00734     }
00735 
00736   // All parameters have been reset.
00737   // We didn't cache the original rhs or matrix for memory reasons,
00738   // but we can restore them to a state consistent solution -
00739   // principle of least surprise.
00740   this->assembly(true, true);
00741   this->rhs->close();
00742   this->matrix->close();
00743   this->assemble_qoi(qoi_indices);
00744 }

std::pair< unsigned int, Real > ImplicitSystem::adjoint_solve ( const QoISet & qoi_indices = QoISet() ) [virtual, inherited]

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

Leave qoi_indices empty to solve all adjoint problems.

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

Reimplemented from System.

Definition at line 342 of file implicit_system.C.

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

00343 {
00344   // Log how long the linear solve takes.
00345   START_LOG("adjoint_solve()", "ImplicitSystem");
00346 
00347   // The forward system should now already be solved.
00348   // Now assemble it's adjoint.
00349 
00350   if (this->assemble_before_solve)
00351     {
00352       // Build the Jacobian
00353       this->assembly(false, true);
00354       this->matrix->close();
00355 
00356       // Take the discrete adjoint
00357       matrix->get_transpose(*matrix);
00358 
00359       // Reset and build the RHS from the QOI derivative
00360       this->assemble_qoi_derivative(qoi_indices);
00361     }
00362 
00363   // The adjoint problem is linear
00364   LinearSolver<Number> *linear_solver = this->get_linear_solver();
00365 
00366   // Our iteration counts and residuals will be sums of the individual
00367   // results
00368   std::pair<unsigned int, Real> solver_params =
00369     this->get_linear_solve_parameters();
00370   std::pair<unsigned int, Real> totalrval = std::make_pair(0,0.0);
00371 
00372   for (unsigned int i=0; i != this->qoi.size(); ++i)
00373     if (qoi_indices.has_index(i))
00374       {
00375         const std::pair<unsigned int, Real> rval =
00376           linear_solver->solve (*matrix, this->add_adjoint_solution(i),
00377                                  this->get_adjoint_rhs(i),
00378                                  solver_params.second,
00379                                  solver_params.first);
00380 
00381         totalrval.first  += rval.first;
00382         totalrval.second += rval.second;
00383       }
00384 
00385   this->release_linear_solver(linear_solver);
00386 
00387   // The linear solver may not have fit our constraints exactly
00388 #ifdef LIBMESH_ENABLE_AMR
00389   for (unsigned int i=0; i != this->qoi.size(); ++i)
00390     if (qoi_indices.has_index(i))
00391       this->get_dof_map().enforce_constraints_exactly
00392         (*this, &this->get_adjoint_solution(i));
00393 #endif
00394 
00395   // Stop logging the nonlinear solve
00396   STOP_LOG("adjoint_solve()", "ImplicitSystem");
00397 
00398   return totalrval;
00399 }

void DifferentiableSystem::assemble ( ) [virtual, inherited]

Prepares matrix and rhs for matrix assembly. Users should not reimplement this

Reimplemented from ImplicitSystem.

Definition at line 102 of file diff_system.C.

References DifferentiableSystem::assembly().

00103 {
00104   this->assembly(true, true);
00105 }

void FEMSystem::assemble_qoi ( const QoISet & indices = QoISet() ) [virtual]

Runs a qoi assembly loop over all elements, and if assemble_qoi_sides is true over all sides.

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

Reimplemented from ExplicitSystem.

Definition at line 486 of file fem_system.C.

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

00487 {
00488   START_LOG("assemble_qoi()", "FEMSystem");
00489 
00490   const MeshBase& mesh = this->get_mesh();
00491 
00492   this->update();
00493 
00494   AutoPtr<DiffContext> con = this->build_context();
00495   FEMContext &_femcontext = libmesh_cast_ref<FEMContext&>(*con);
00496 
00497   // the quantity of interest is assumed to be a sum of element and
00498   // side terms
00499   for (unsigned int i=0; i != qoi.size(); ++i)
00500     if (qoi_indices.has_index(i))
00501       qoi[i] = 0;
00502 
00503   // Loop over every active mesh element on this processor
00504   MeshBase::const_element_iterator el =
00505     mesh.active_local_elements_begin();
00506   const MeshBase::const_element_iterator end_el =
00507     mesh.active_local_elements_end();
00508 
00509   for ( ; el != end_el; ++el)
00510     {
00511       _femcontext.reinit(*this, *el);
00512 
00513       this->element_qoi(_femcontext, qoi_indices);
00514 
00515       for (_femcontext.side = 0;
00516            _femcontext.side != _femcontext.elem->n_sides();
00517            ++_femcontext.side)
00518         {
00519           // Don't compute on non-boundary sides unless requested
00520           if (!assemble_qoi_sides ||
00521               (!compute_internal_sides &&
00522                _femcontext.elem->neighbor(_femcontext.side) != NULL))
00523             continue;
00524 
00525           std::map<FEType, FEBase *>::iterator fe_end =
00526             _femcontext.side_fe.end();
00527           for (std::map<FEType, FEBase *>::iterator i =
00528             _femcontext.side_fe.begin();
00529                i != fe_end; ++i)
00530             {
00531               i->second->reinit(_femcontext.elem, _femcontext.side);
00532             }
00533 
00534           this->side_qoi(_femcontext, qoi_indices);
00535         }
00536     }
00537 
00538   Parallel::sum(_femcontext.elem_qoi);
00539 
00540   this->qoi = _femcontext.elem_qoi;
00541 
00542   STOP_LOG("assemble_qoi()", "FEMSystem");
00543 }

void FEMSystem::assemble_qoi_derivative ( const QoISet & indices = QoISet() ) [virtual]

Runs a qoi derivative assembly loop over all elements, and if assemble_qoi_sides is true over all sides.

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

Reimplemented from ExplicitSystem.

Definition at line 547 of file fem_system.C.

References MeshBase::active_local_elements_begin(), MeshBase::active_local_elements_end(), System::add_adjoint_rhs(), DifferentiableSystem::assemble_qoi_sides, build_context(), DifferentiableSystem::compute_internal_sides, DofMap::constrain_element_vector(), DiffContext::dof_indices, FEMContext::elem, DiffContext::elem_qoi_derivative, DifferentiableSystem::element_qoi_derivative(), System::get_adjoint_rhs(), System::get_dof_map(), System::get_mesh(), QoISet::has_index(), mesh, Elem::n_sides(), Elem::neighbor(), System::qoi, FEMContext::reinit(), FEMContext::side, FEMContext::side_fe, DifferentiableSystem::side_qoi_derivative(), and System::update().

00548 {
00549   START_LOG("assemble_qoi_derivative()", "FEMSystem");
00550 
00551   const MeshBase& mesh = this->get_mesh();
00552 
00553   this->update();
00554 
00555   AutoPtr<DiffContext> con = this->build_context();
00556   FEMContext &_femcontext = libmesh_cast_ref<FEMContext&>(*con);
00557 
00558   // The quantity of interest derivative assembly accumulates on
00559   // initially zero vectors
00560   for (unsigned int i=0; i != qoi.size(); ++i)
00561     if (qoi_indices.has_index(i))
00562       this->add_adjoint_rhs(i).zero();
00563 
00564   // Loop over every active mesh element on this processor
00565   MeshBase::const_element_iterator el =
00566     mesh.active_local_elements_begin();
00567   const MeshBase::const_element_iterator end_el =
00568     mesh.active_local_elements_end();
00569 
00570   for ( ; el != end_el; ++el)
00571     {
00572       _femcontext.reinit(*this, *el);
00573 
00574       this->element_qoi_derivative(_femcontext, qoi_indices);
00575 
00576       for (_femcontext.side = 0;
00577            _femcontext.side != _femcontext.elem->n_sides();
00578            ++_femcontext.side)
00579         {
00580           // Don't compute on non-boundary sides unless requested
00581           if (!assemble_qoi_sides ||
00582               (!compute_internal_sides &&
00583                _femcontext.elem->neighbor(_femcontext.side) != NULL))
00584             continue;
00585 
00586           std::map<FEType, FEBase *>::iterator fe_end =
00587             _femcontext.side_fe.end();
00588           for (std::map<FEType, FEBase *>::iterator i =
00589             _femcontext.side_fe.begin();
00590                i != fe_end; ++i)
00591             {
00592               i->second->reinit(_femcontext.elem, _femcontext.side);
00593             }
00594 
00595           this->side_qoi_derivative(_femcontext, qoi_indices);
00596         }
00597 
00598       // We need some unmodified indices to use for constraining
00599       // multiple vector
00600       // FIXME - there should be a DofMap::constrain_element_vectors
00601       // to do this more efficiently
00602       std::vector<unsigned int> original_dofs = _femcontext.dof_indices;
00603 
00604       for (unsigned int i=0; i != qoi.size(); ++i)
00605         if (qoi_indices.has_index(i))
00606           {
00607             _femcontext.dof_indices = original_dofs;
00608             this->get_dof_map().constrain_element_vector
00609               (_femcontext.elem_qoi_derivative[i], _femcontext.dof_indices, false);
00610 
00611             this->get_adjoint_rhs(i).add_vector
00612               (_femcontext.elem_qoi_derivative[i], _femcontext.dof_indices);
00613           }
00614     }
00615 
00616   STOP_LOG("assemble_qoi_derivative()", "FEMSystem");
00617 }

void ImplicitSystem::assemble_residual_derivatives ( const ParameterVector & parameters ) [virtual, inherited]

Residual parameter derivative function.

Uses finite differences by default.

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

Can be overloaded in derived classes.

Reimplemented from System.

Definition at line 622 of file implicit_system.C.

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

Referenced by ImplicitSystem::sensitivity_solve().

00623 {
00624   const unsigned int Np = parameters.size();
00625   Real deltap = TOLERANCE;
00626 
00627   for (unsigned int p=0; p != Np; ++p)
00628     {
00629       NumericVector<Number> &sensitivity_rhs = this->add_sensitivity_rhs(p);
00630 
00631       // Approximate -(partial R / partial p) by
00632       // (R(p-dp) - R(p+dp)) / (2*dp)
00633 
00634       Number old_parameter = *parameters[p];
00635       *parameters[p] -= deltap;
00636 
00637       this->assembly(true, false);
00638       this->rhs->close();
00639       sensitivity_rhs = *this->rhs;
00640 
00641       *parameters[p] = old_parameter + deltap;
00642 
00643       this->assembly(true, false);
00644       this->rhs->close();
00645 
00646       sensitivity_rhs -= *this->rhs;
00647       sensitivity_rhs /= (2*deltap);
00648       sensitivity_rhs.close();
00649 
00650       *parameters[p] = old_parameter;
00651     }
00652 }

void FEMSystem::assembly	(	bool	get_residual,
		bool	get_jacobian
	)			`[virtual]`

Prepares matrix or rhs for matrix assembly. Users may reimplement this to add pre- or post-assembly code before or after calling FEMSystem::assembly()

Implements DifferentiableSystem.

Definition at line 78 of file fem_system.C.

Referenced by ContinuationSystem::continuation_solve(), and ContinuationSystem::solve_tangent().

00079 {
00080   libmesh_assert(get_residual || get_jacobian);
00081   std::string log_name;
00082   if (get_residual && get_jacobian)
00083     log_name = "assembly()";
00084   else if (get_residual)
00085     log_name = "assembly(get_residual)";
00086   else
00087     log_name = "assembly(get_jacobian)";
00088 
00089   START_LOG(log_name, "FEMSystem");
00090 
00091   const MeshBase& mesh = this->get_mesh();
00092 
00093 //  this->get_vector("_nonlinear_solution").localize
00094 //    (*current_local_nonlinear_solution,
00095 //     dof_map.get_send_list());
00096   this->update();
00097 
00098   if (print_solution_norms)
00099     {
00100 //      this->get_vector("_nonlinear_solution").close();
00101       this->solution->close();
00102       std::cout << "|U| = "
00103 //                << this->get_vector("_nonlinear_solution").l1_norm()
00104                 << this->solution->l1_norm()
00105                 << std::endl;
00106     }
00107   if (print_solutions)
00108     {
00109       unsigned int old_precision = std::cout.precision();
00110       std::cout.precision(16);
00111 //      std::cout << "U = [" << this->get_vector("_nonlinear_solution")
00112       std::cout << "U = [" << *(this->solution)
00113                 << "];" << std::endl;
00114       std::cout.precision(old_precision);
00115     }
00116 
00117   // Is this definitely necessary? [RHS]
00118   if (get_jacobian)
00119     matrix->zero();
00120   if (get_residual)
00121     rhs->zero();
00122 
00123   // Stupid C++ lets you set *Real* verify_analytic_jacobians = true!
00124   if (verify_analytic_jacobians > 0.5)
00125     {
00126       std::cerr << "WARNING!  verify_analytic_jacobians was set "
00127                 << "to absurdly large value of "
00128                 << verify_analytic_jacobians << std::endl;
00129       std::cerr << "Resetting to 1e-6!" << std::endl;
00130       verify_analytic_jacobians = 1e-6;
00131     }
00132 
00133   // In time-dependent problems, the nonlinear function we're trying
00134   // to solve at each timestep may depend on the particular solver
00135   // we're using
00136   libmesh_assert (time_solver.get() != NULL);
00137 
00138   AutoPtr<DiffContext> con = this->build_context();
00139   FEMContext &_femcontext = libmesh_cast_ref<FEMContext&>(*con);
00140 
00141   // Build the residual and jacobian contributions on every active
00142   // mesh element on this processor
00143   MeshBase::const_element_iterator el =
00144     mesh.active_local_elements_begin();
00145   const MeshBase::const_element_iterator end_el =
00146     mesh.active_local_elements_end();
00147 
00148   for ( ; el != end_el; ++el)
00149     {
00150       _femcontext.reinit(*this, *el);
00151 
00152       bool jacobian_computed =
00153         time_solver->element_residual(get_jacobian, _femcontext);
00154 
00155       // Compute a numeric jacobian if we have to
00156       if (get_jacobian && !jacobian_computed)
00157         {
00158           // Make sure we didn't compute a jacobian and lie about it
00159           libmesh_assert(_femcontext.elem_jacobian.l1_norm() == 0.0);
00160           // Logging of numerical jacobians is done separately
00161           this->numerical_elem_jacobian(_femcontext);
00162         }
00163 
00164       // Compute a numeric jacobian if we're asked to verify the
00165       // analytic jacobian we got
00166       if (get_jacobian && jacobian_computed &&
00167           verify_analytic_jacobians != 0.0)
00168         {
00169           DenseMatrix<Number> analytic_jacobian(_femcontext.elem_jacobian);
00170 
00171           _femcontext.elem_jacobian.zero();
00172           // Logging of numerical jacobians is done separately
00173           this->numerical_elem_jacobian(_femcontext);
00174 
00175           Real analytic_norm = analytic_jacobian.l1_norm();
00176           Real numerical_norm = _femcontext.elem_jacobian.l1_norm();
00177 
00178           // If we can continue, we'll probably prefer the analytic jacobian
00179           analytic_jacobian.swap(_femcontext.elem_jacobian);
00180 
00181           // The matrix "analytic_jacobian" will now hold the error matrix
00182           analytic_jacobian.add(-1.0, _femcontext.elem_jacobian);
00183           Real error_norm = analytic_jacobian.l1_norm();
00184 
00185           Real relative_error = error_norm /
00186                                 std::max(analytic_norm, numerical_norm);
00187 
00188           if (relative_error > verify_analytic_jacobians)
00189             {
00190               std::cerr << "Relative error " << relative_error
00191                         << " detected in analytic jacobian on element "
00192                         << _femcontext.elem->id() << '!' << std::endl;
00193 
00194               unsigned int old_precision = std::cout.precision();
00195               std::cout.precision(16);
00196               std::cout << "J_analytic " << _femcontext.elem->id() << " = "
00197                         << _femcontext.elem_jacobian << std::endl;
00198               analytic_jacobian.add(1.0, _femcontext.elem_jacobian);
00199               std::cout << "J_numeric " << _femcontext.elem->id() << " = "
00200                         << analytic_jacobian << std::endl;
00201 
00202               std::cout.precision(old_precision);
00203 
00204               libmesh_error();
00205             }
00206         }
00207 
00208       for (_femcontext.side = 0;
00209            _femcontext.side != _femcontext.elem->n_sides();
00210            ++_femcontext.side)
00211         {
00212           // Don't compute on non-boundary sides unless requested
00213           if (!compute_internal_sides && 
00214               _femcontext.elem->neighbor(_femcontext.side) != NULL)
00215             continue;
00216 
00217           // Any mesh movement has already been done (and restored,
00218           // if the TimeSolver isn't broken), but
00219           // reinitializing the side FE objects is still necessary
00220           _femcontext.elem_side_fe_reinit();
00221 
00222           DenseMatrix<Number> old_jacobian;
00223           // If we're in DEBUG mode, we should always verify that the
00224           // user's side_residual function doesn't alter our existing
00225           // jacobian and then lie about it
00226 #ifndef DEBUG
00227           // Even if we're not in DEBUG mode, when we're verifying
00228           // analytic jacobians we'll want to verify each side's
00229           // jacobian contribution separately
00230           if (verify_analytic_jacobians != 0.0 && get_jacobian)
00231 #endif // ifndef DEBUG
00232             {
00233               old_jacobian = _femcontext.elem_jacobian;
00234               _femcontext.elem_jacobian.zero();
00235             }
00236           jacobian_computed =
00237             time_solver->side_residual(get_jacobian, _femcontext);
00238 
00239           // Compute a numeric jacobian if we have to
00240           if (get_jacobian && !jacobian_computed)
00241             {
00242               // In DEBUG mode, we've already set elem_jacobian == 0,
00243               // so we can make sure side_residual didn't compute a
00244               // jacobian and lie about it
00245 #ifdef DEBUG
00246               libmesh_assert(_femcontext.elem_jacobian.l1_norm() == 0.0);
00247 #endif
00248               // Logging of numerical jacobians is done separately
00249               this->numerical_side_jacobian(_femcontext);
00250 
00251               // If we're in DEBUG mode or if
00252               // verify_analytic_jacobians is on, we've moved
00253               // elem_jacobian's accumulated values into old_jacobian.
00254               // Now let's add them back.
00255 #ifndef DEBUG
00256               if (verify_analytic_jacobians != 0.0)
00257 #endif // ifndef DEBUG
00258                 _femcontext.elem_jacobian += old_jacobian;
00259             }
00260 
00261           // Compute a numeric jacobian if we're asked to verify the
00262           // analytic jacobian we got
00263           if (get_jacobian && jacobian_computed &&
00264               verify_analytic_jacobians != 0.0)
00265             {
00266               DenseMatrix<Number> analytic_jacobian(_femcontext.elem_jacobian);
00267 
00268               _femcontext.elem_jacobian.zero();
00269               // Logging of numerical jacobians is done separately
00270               this->numerical_side_jacobian(_femcontext);
00271 
00272               Real analytic_norm = analytic_jacobian.l1_norm();
00273               Real numerical_norm = _femcontext.elem_jacobian.l1_norm();
00274 
00275               // If we can continue, we'll probably prefer the analytic jacobian
00276               analytic_jacobian.swap(_femcontext.elem_jacobian);
00277 
00278               // The matrix "analytic_jacobian" will now hold the error matrix
00279               analytic_jacobian.add(-1.0, _femcontext.elem_jacobian);
00280               Real error_norm = analytic_jacobian.l1_norm();
00281 
00282               Real relative_error = error_norm /
00283                                     std::max(analytic_norm, numerical_norm);
00284 
00285               if (relative_error > verify_analytic_jacobians)
00286                 {
00287                   std::cerr << "Relative error " << relative_error
00288                             << " detected in analytic jacobian on element "
00289                             << _femcontext.elem->id()
00290                             << ", side "
00291                             << _femcontext.side << '!' << std::endl;
00292 
00293                   unsigned int old_precision = std::cout.precision();
00294                   std::cout.precision(16);
00295                   std::cout << "J_analytic " << _femcontext.elem->id() << " = "
00296                             << _femcontext.elem_jacobian << std::endl;
00297                   analytic_jacobian.add(1.0, _femcontext.elem_jacobian);
00298                   std::cout << "J_numeric " << _femcontext.elem->id() << " = "
00299                             << analytic_jacobian << std::endl;
00300                   std::cout.precision(old_precision);
00301 
00302                   libmesh_error();
00303                 }
00304               // Once we've verified a side, we'll want to add back the
00305               // rest of the accumulated jacobian
00306               _femcontext.elem_jacobian += old_jacobian;
00307             }
00308           // In DEBUG mode, we've set elem_jacobian == 0, and we
00309           // may still need to add the old jacobian back
00310 #ifdef DEBUG
00311           if (get_jacobian && jacobian_computed &&
00312               verify_analytic_jacobians == 0.0)
00313             {
00314               _femcontext.elem_jacobian += old_jacobian;
00315             }
00316 #endif // ifdef DEBUG
00317         }
00318 
00319 #ifdef LIBMESH_ENABLE_AMR
00320       // We turn off the asymmetric constraint application;
00321       // enforce_constraints_exactly() should be called in the solver
00322       if (get_residual && get_jacobian)
00323         this->get_dof_map().constrain_element_matrix_and_vector
00324           (_femcontext.elem_jacobian, _femcontext.elem_residual,
00325            _femcontext.dof_indices, false);
00326       else if (get_residual)
00327         this->get_dof_map().constrain_element_vector
00328           (_femcontext.elem_residual, _femcontext.dof_indices, false);
00329       else if (get_jacobian)
00330         this->get_dof_map().constrain_element_matrix
00331           (_femcontext.elem_jacobian, _femcontext.dof_indices, false);
00332 #endif // #ifdef LIBMESH_ENABLE_AMR
00333 
00334       if (get_jacobian && print_element_jacobians)
00335         {
00336           unsigned int old_precision = std::cout.precision();
00337           std::cout.precision(16);
00338           std::cout << "J_elem " << _femcontext.elem->id() << " = "
00339                     << _femcontext.elem_jacobian << std::endl;
00340           std::cout.precision(old_precision);
00341         }
00342 
00343       if (get_jacobian)
00344         this->matrix->add_matrix (_femcontext.elem_jacobian,
00345                                   _femcontext.dof_indices);
00346       if (get_residual)
00347         this->rhs->add_vector (_femcontext.elem_residual,
00348                                _femcontext.dof_indices);
00349     }
00350 
00351 
00352   if (get_residual && print_residual_norms)
00353     {
00354       this->rhs->close();
00355       std::cout << "|F| = " << this->rhs->l1_norm() << std::endl;
00356     }
00357   if (get_residual && print_residuals)
00358     {
00359       this->rhs->close();
00360       unsigned int old_precision = std::cout.precision();
00361       std::cout.precision(16);
00362       std::cout << "F = [" << *(this->rhs) << "];" << std::endl;
00363       std::cout.precision(old_precision);
00364     }
00365   if (get_jacobian && print_jacobian_norms)
00366     {
00367       this->matrix->close();
00368       std::cout << "|J| = " << this->matrix->l1_norm() << std::endl;
00369     }
00370   if (get_jacobian && print_jacobians)
00371     {
00372       this->matrix->close();
00373       unsigned int old_precision = std::cout.precision();
00374       std::cout.precision(16);
00375       std::cout << "J = [" << *(this->matrix) << "];" << std::endl;
00376       std::cout.precision(old_precision);
00377     }
00378   STOP_LOG(log_name, "FEMSystem");
00379 }

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

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

Definition at line 1316 of file system.C.

References System::_assemble_system.

01318 {
01319   libmesh_assert (fptr != NULL);
01320   
01321   _assemble_system = fptr;  
01322 }

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

Register a user function for imposing constraints.

Definition at line 1326 of file system.C.

References System::_constrain_system.

01328 {
01329   libmesh_assert (fptr != NULL);
01330   
01331   _constrain_system = fptr;  
01332 }

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

Register a user function to use in initializing the system.

Definition at line 1306 of file system.C.

References System::_init_system.

01308 {
01309   libmesh_assert (fptr != NULL);
01310   
01311   _init_system = fptr;
01312 }

void System::attach_QOI_derivative ( void fptrEquationSystems &es, const std::string &name, const QoISet &qoi_indices ) [inherited]

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

void System::attach_QOI_function ( void fptrEquationSystems &es, const std::string &name, const QoISet &qoi_indices ) [inherited]

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

AutoPtr< DiffContext > FEMSystem::build_context ( ) [virtual]

Builds a FEMContext object with enough information to do evaluations on each element.

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

Reimplemented from DifferentiableSystem.

Definition at line 742 of file fem_system.C.

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

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

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

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

Definition at line 1077 of file system.C.

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

01079 {
01080   // This function must be run on all processors at once
01081   parallel_only();
01082 
01083   START_LOG ("calculate_norm()", "System");
01084 
01085   // Zero the norm before summation
01086   Real v_norm = 0.;
01087 
01088   if (norm.is_discrete())
01089     {
01090       STOP_LOG ("calculate_norm()", "System");
01091       //Check to see if all weights are 1.0
01092       unsigned int check_var = 0;
01093       for (; check_var != this->n_vars(); ++check_var)
01094         if(norm.weight(check_var) != 1.0)
01095           break;
01096 
01097       //All weights were 1.0 so just do the full vector discrete norm
01098       if(check_var == this->n_vars())
01099         {
01100           FEMNormType norm_type = norm.type(0);
01101           
01102           if(norm_type == DISCRETE_L1)
01103             return v.l1_norm();
01104           if(norm_type == DISCRETE_L2)
01105             return v.l2_norm();
01106           if(norm_type == DISCRETE_L_INF)
01107             return v.linfty_norm();
01108           else
01109             libmesh_error();
01110         }
01111 
01112       for (unsigned int var=0; var != this->n_vars(); ++var)
01113         {
01114           // Skip any variables we don't need to integrate
01115           if (norm.weight(var) == 0.0)
01116             continue;
01117 
01118           v_norm += norm.weight(var) * discrete_var_norm(v, var, norm.type(var));
01119         }
01120 
01121       return v_norm;
01122     }
01123 
01124   // Localize the potentially parallel vector
01125   AutoPtr<NumericVector<Number> > local_v = NumericVector<Number>::build();
01126   local_v->init(v.size(), true, SERIAL);
01127   v.localize (*local_v, _dof_map->get_send_list()); 
01128 
01129   unsigned int dim = this->get_mesh().mesh_dimension();
01130 
01131   // Loop over all variables
01132   for (unsigned int var=0; var != this->n_vars(); ++var)
01133     {
01134       // Skip any variables we don't need to integrate
01135       if (norm.weight(var) == 0.0)
01136         continue;
01137 
01138       const FEType& fe_type = this->get_dof_map().variable_type(var);
01139       AutoPtr<QBase> qrule =
01140         fe_type.default_quadrature_rule (dim);
01141       AutoPtr<FEBase> fe
01142         (FEBase::build(dim, fe_type));
01143       fe->attach_quadrature_rule (qrule.get());
01144 
01145       const std::vector<Real>&               JxW = fe->get_JxW();
01146       const std::vector<std::vector<Real> >* phi = NULL;
01147       if (norm.type(var) == H1 ||
01148           norm.type(var) == H2 ||
01149           norm.type(var) == L2)
01150         phi = &(fe->get_phi());
01151 
01152       const std::vector<std::vector<RealGradient> >* dphi = NULL;
01153       if (norm.type(var) == H1 ||
01154           norm.type(var) == H2 ||
01155           norm.type(var) == H1_SEMINORM)
01156         dphi = &(fe->get_dphi());
01157 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
01158       const std::vector<std::vector<RealTensor> >*   d2phi = NULL;
01159       if (norm.type(var) == H2 ||
01160           norm.type(var) == H2_SEMINORM)
01161         d2phi = &(fe->get_d2phi());
01162 #endif
01163 
01164       std::vector<unsigned int> dof_indices;
01165 
01166       // Begin the loop over the elements
01167       MeshBase::const_element_iterator       el     =
01168         this->get_mesh().active_local_elements_begin();
01169       const MeshBase::const_element_iterator end_el =
01170         this->get_mesh().active_local_elements_end();
01171 
01172       for ( ; el != end_el; ++el)
01173         {
01174           const Elem* elem = *el;
01175 
01176           fe->reinit (elem);
01177 
01178           this->get_dof_map().dof_indices (elem, dof_indices, var);
01179 
01180           const unsigned int n_qp = qrule->n_points();
01181 
01182           const unsigned int n_sf = dof_indices.size();
01183 
01184           // Begin the loop over the Quadrature points.
01185           for (unsigned int qp=0; qp<n_qp; qp++)
01186             {
01187               if (norm.type(var) == H1 ||
01188                   norm.type(var) == H2 ||
01189                   norm.type(var) == L2)
01190                 {
01191                   Number u_h = 0.;
01192                   for (unsigned int i=0; i != n_sf; ++i)
01193                     u_h += (*phi)[i][qp] * (*local_v)(dof_indices[i]);
01194                   v_norm += norm.weight_sq(var) *
01195                             JxW[qp] * libmesh_norm(u_h);
01196                 }
01197 
01198               if (norm.type(var) == H1 ||
01199                   norm.type(var) == H2 ||
01200                   norm.type(var) == H1_SEMINORM)
01201                 {
01202                   Gradient grad_u_h;
01203                   for (unsigned int i=0; i != n_sf; ++i)
01204                     grad_u_h.add_scaled((*dphi)[i][qp], (*local_v)(dof_indices[i]));
01205                   v_norm += norm.weight_sq(var) *
01206                             JxW[qp] * grad_u_h.size_sq();
01207                 }
01208 
01209 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
01210               if (norm.type(var) == H2 ||
01211                   norm.type(var) == H2_SEMINORM)
01212                 {
01213                   Tensor hess_u_h;
01214                   for (unsigned int i=0; i != n_sf; ++i)
01215                     hess_u_h.add_scaled((*d2phi)[i][qp], (*local_v)(dof_indices[i]));
01216                   v_norm += norm.weight_sq(var) *
01217                             JxW[qp] * hess_u_h.size_sq();
01218                 }
01219 #endif
01220             }
01221         }
01222     }
01223 
01224   Parallel::sum(v_norm);
01225 
01226   STOP_LOG ("calculate_norm()", "System");
01227 
01228   return std::sqrt(v_norm);
01229 }

Real System::calculate_norm	(	NumericVector< Number > &	v,
		unsigned int	var = `0`,
		FEMNormType	norm_type = `L2`
	)			const `[inherited]`

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

Definition at line 1056 of file system.C.

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

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

01059 {
01060   //short circuit to save time
01061   if(norm_type == DISCRETE_L1 ||
01062      norm_type == DISCRETE_L2 ||
01063      norm_type == DISCRETE_L_INF)
01064     return discrete_var_norm(v,var,norm_type);
01065 
01066   // Not a discrete norm
01067   std::vector<FEMNormType> norms(this->n_vars(), L2);
01068   std::vector<Real> weights(this->n_vars(), 0.0);
01069   norms[var] = norm_type;
01070   weights[var] = 1.0;
01071   Real val = this->calculate_norm(v, SystemNorm(norms, weights));
01072   return val;
01073 }

void FEMSystem::clear ( ) [virtual]

Clear all the data structures associated with the system.

Reimplemented from DifferentiableSystem.

Reimplemented in ContinuationSystem.

Definition at line 64 of file fem_system.C.

References DifferentiableSystem::clear().

Referenced by ContinuationSystem::clear(), and ~FEMSystem().

00065 {
00066   Parent::clear();
00067 }

bool System::compare	(	const System &	other_system,
		const Real	threshold,
		const bool	verbose
	)			const `[virtual, inherited]`

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

Definition at line 399 of file system.C.

References System::_can_add_vectors, System::_sys_name, System::_vectors, System::get_vector(), System::n_vectors(), System::name(), and System::solution.

Referenced by EquationSystems::compare().

00402 {
00403   // we do not care for matrices, but for vectors
00404   libmesh_assert (!_can_add_vectors);
00405   libmesh_assert (!other_system._can_add_vectors);
00406 
00407   if (verbose)
00408     {
00409       std::cout << "  Systems \"" << _sys_name << "\"" << std::endl;
00410       std::cout << "   comparing matrices not supported." << std::endl;
00411       std::cout << "   comparing names...";
00412     }
00413 
00414   // compare the name: 0 means identical
00415   const int name_result = _sys_name.compare(other_system.name());
00416   if (verbose)
00417     {
00418       if (name_result == 0)
00419         std::cout << " identical." << std::endl;
00420       else
00421         std::cout << "  names not identical." << std::endl;
00422       std::cout << "   comparing solution vector...";
00423     }
00424 
00425 
00426   // compare the solution: -1 means identical
00427   const int solu_result = solution->compare (*other_system.solution.get(),
00428                                              threshold);
00429 
00430   if (verbose)
00431     {
00432       if (solu_result == -1)
00433         std::cout << " identical up to threshold." << std::endl;
00434       else
00435         std::cout << "  first difference occured at index = " 
00436                   << solu_result << "." << std::endl;
00437     }
00438 
00439 
00440   // safety check, whether we handle at least the same number
00441   // of vectors
00442   std::vector<int> ov_result;
00443 
00444   if (this->n_vectors() != other_system.n_vectors())
00445     {
00446       if (verbose)
00447         {
00448           std::cout << "   Fatal difference. This system handles " 
00449                     << this->n_vectors() << " add'l vectors," << std::endl
00450                     << "   while the other system handles "
00451                     << other_system.n_vectors() 
00452                     << " add'l vectors." << std::endl
00453                     << "   Aborting comparison." << std::endl;
00454         }
00455       return false;
00456     }
00457   else if (this->n_vectors() == 0)
00458     {
00459       // there are no additional vectors...
00460       ov_result.clear ();
00461     }
00462   else
00463     {
00464       // compare other vectors
00465       for (const_vectors_iterator pos = _vectors.begin();
00466            pos != _vectors.end(); ++pos)
00467         {
00468           if (verbose)
00469               std::cout << "   comparing vector \""
00470                         << pos->first << "\" ...";
00471 
00472           // assume they have the same name
00473           const NumericVector<Number>& other_system_vector = 
00474               other_system.get_vector(pos->first);
00475 
00476           ov_result.push_back(pos->second->compare (other_system_vector,
00477                                                     threshold));
00478 
00479           if (verbose)
00480             {
00481               if (ov_result[ov_result.size()-1] == -1)
00482                 std::cout << " identical up to threshold." << std::endl;
00483               else
00484                 std::cout << " first difference occured at" << std::endl
00485                           << "   index = " << ov_result[ov_result.size()-1] << "." << std::endl;
00486             }
00487 
00488         }
00489 
00490     } // finished comparing additional vectors
00491 
00492 
00493   bool overall_result;
00494        
00495   // sum up the results
00496   if ((name_result==0) && (solu_result==-1))
00497     {
00498       if (ov_result.size()==0)
00499         overall_result = true;
00500       else
00501         {
00502           bool ov_identical;
00503           unsigned int n    = 0;
00504           do
00505             {
00506               ov_identical = (ov_result[n]==-1);
00507               n++;
00508             }
00509           while (ov_identical && n<ov_result.size());
00510           overall_result = ov_identical;
00511         }
00512     }
00513   else
00514     overall_result = false;
00515 
00516   if (verbose)
00517     {
00518       std::cout << "   finished comparisons, ";
00519       if (overall_result)
00520         std::cout << "found no differences." << std::endl << std::endl;
00521       else 
00522         std::cout << "found differences." << std::endl << std::endl;
00523     }
00524           
00525   return overall_result;
00526 }

Number System::current_solution ( const unsigned int global_dof_number ) const [inherited]

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

Definition at line 117 of file system.C.

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

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

00118 {
00119   // Check the sizes
00120   libmesh_assert (global_dof_number < _dof_map->n_dofs());
00121   libmesh_assert (global_dof_number < current_local_solution->size());
00122    
00123   return (*current_local_solution)(global_dof_number);
00124 }

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

Deactivates the system. Only active systems are solved.

Definition at line 1388 of file system.h.

References System::_active.

01389 {
01390   _active = false;
01391 }

virtual bool DifferentiableSystem::element_constraint	(	bool	request_jacobian,
		DiffContext &
	)			`[inline, virtual, inherited]`

Adds the constraint contribution on elem to elem_residual. If this method receives request_jacobian = true, then it should compute elem_jacobian and return true if possible. If elem_jacobian has not been computed then the method should return false.

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

Definition at line 157 of file diff_system.h.

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

00158                                                   {
00159     return request_jacobian;
00160   }

virtual void DifferentiableSystem::element_postprocess ( DiffContext & ) [inline, virtual, inherited]

Does any work that needs to be done on elem in a postprocessing loop.

Definition at line 225 of file diff_system.h.

Referenced by postprocess().

00225 {}

virtual void DifferentiableSystem::element_qoi	(	DiffContext &	,
		const QoISet &
	)			`[inline, virtual, inherited]`

Does any work that needs to be done on elem in a quantity of interest assembly loop, outputting to elem_qoi.

Only qois included in the supplied QoISet need to be assembled.

Definition at line 240 of file diff_system.h.

Referenced by assemble_qoi().

00242     {}

virtual void DifferentiableSystem::element_qoi_derivative	(	DiffContext &	,
		const QoISet &
	)			`[inline, virtual, inherited]`

Does any work that needs to be done on elem in a quantity of interest derivative assembly loop, outputting to elem_qoi_derivative

Only qois included in the supplied QoISet need their derivatives assembled.

Definition at line 252 of file diff_system.h.

Referenced by assemble_qoi_derivative().

00254     {}

virtual bool DifferentiableSystem::element_time_derivative	(	bool	request_jacobian,
		DiffContext &
	)			`[inline, virtual, inherited]`

Adds the time derivative contribution on elem to elem_residual. If this method receives request_jacobian = true, then it should compute elem_jacobian and return true if possible. If elem_jacobian has not been computed then the method should return false.

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

Definition at line 140 of file diff_system.h.

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

00141                                                        {
00142     return request_jacobian;
00143   }

bool FEMSystem::eulerian_residual	(	bool	request_jacobian,
		DiffContext &	context
	)			`[virtual]`

Adds a pseudo-convection contribution on elem to elem_residual, if the nodes of elem are being translated by a moving mesh.

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

Reimplemented from DifferentiableSystem.

Definition at line 864 of file fem_system.C.

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

00866 {
00867   // Only calculate a mesh movement residual if it's necessary
00868   if (_mesh_sys == libMesh::invalid_uint)
00869     return request_jacobian;
00870 
00871   FEMContext &context = libmesh_cast_ref<FEMContext&>(c);
00872 
00873   // This function only supports fully coupled mesh motion for now
00874   libmesh_assert(_mesh_sys == this->number());
00875 
00876   unsigned int n_qpoints = context.element_qrule->n_points();
00877 
00878   const unsigned int n_x_dofs = (_mesh_x_var == libMesh::invalid_uint) ?
00879                                 0 : context.dof_indices_var[_mesh_x_var].size();
00880   const unsigned int n_y_dofs = (_mesh_y_var == libMesh::invalid_uint) ?
00881                                 0 : context.dof_indices_var[_mesh_y_var].size();
00882   const unsigned int n_z_dofs = (_mesh_z_var == libMesh::invalid_uint) ?
00883                                 0 : context.dof_indices_var[_mesh_z_var].size();
00884 
00885   const unsigned int mesh_xyz_var = n_x_dofs ? _mesh_x_var :
00886                                    (n_y_dofs ? _mesh_y_var :
00887                                    (n_z_dofs ? _mesh_z_var :
00888                                     libMesh::invalid_uint));
00889 
00890   // If we're our own _mesh_sys, we'd better be in charge of
00891   // at least one coordinate, and we'd better have the same
00892   // FE type for all coordinates we are in charge of
00893   libmesh_assert(mesh_xyz_var != libMesh::invalid_uint);
00894   libmesh_assert(!n_x_dofs || context.element_fe_var[_mesh_x_var] ==
00895                               context.element_fe_var[mesh_xyz_var]);
00896   libmesh_assert(!n_y_dofs || context.element_fe_var[_mesh_y_var] ==
00897                               context.element_fe_var[mesh_xyz_var]);
00898   libmesh_assert(!n_z_dofs || context.element_fe_var[_mesh_z_var] ==
00899                               context.element_fe_var[mesh_xyz_var]);
00900 
00901   const std::vector<std::vector<Real> >     &psi =
00902     context.element_fe_var[mesh_xyz_var]->get_phi();
00903 
00904   for (unsigned int var = 0; var != this->n_vars(); ++var)
00905     {
00906       // Mesh motion only affects time-evolving variables
00907       if (!_time_evolving[var])
00908         continue;
00909 
00910       // The mesh coordinate variables themselves are Lagrangian,
00911       // not Eulerian, and no convective term is desired.
00912       if (_mesh_sys == this->number() &&
00913           (var == _mesh_x_var ||
00914            var == _mesh_y_var ||
00915            var == _mesh_z_var))
00916         continue;
00917 
00918       // Some of this code currently relies on the assumption that
00919       // we can pull mesh coordinate data from our own system
00920       if (_mesh_sys != this->number())
00921         libmesh_not_implemented();
00922 
00923       // This residual should only be called by unsteady solvers:
00924       // if the mesh is steady, there's no mesh convection term!
00925       UnsteadySolver *unsteady =
00926         dynamic_cast<UnsteadySolver *>(this->time_solver.get());
00927       if (!unsteady)
00928         return request_jacobian;
00929 
00930       const std::vector<Real> &JxW = 
00931         context.element_fe_var[var]->get_JxW();
00932 
00933       const std::vector<std::vector<Real> >     &phi =
00934         context.element_fe_var[var]->get_phi();
00935 
00936       const std::vector<std::vector<RealGradient> > &dphi =
00937         context.element_fe_var[var]->get_dphi();
00938 
00939       const unsigned int n_u_dofs = context.dof_indices_var[var].size();
00940 
00941       DenseSubVector<Number> &Fu = *context.elem_subresiduals[var];
00942       DenseSubMatrix<Number> &Kuu = *context.elem_subjacobians[var][var];
00943 
00944       DenseSubMatrix<Number> *Kux = n_x_dofs ?
00945         context.elem_subjacobians[var][_mesh_x_var] : NULL;
00946       DenseSubMatrix<Number> *Kuy = n_y_dofs ?
00947         context.elem_subjacobians[var][_mesh_y_var] : NULL;
00948       DenseSubMatrix<Number> *Kuz = n_z_dofs ?
00949         context.elem_subjacobians[var][_mesh_z_var] : NULL;
00950 
00951       std::vector<Real> delta_x(n_x_dofs, 0.);
00952       std::vector<Real> delta_y(n_y_dofs, 0.);
00953       std::vector<Real> delta_z(n_z_dofs, 0.);
00954 
00955       for (unsigned int i = 0; i != n_x_dofs; ++i)
00956         {
00957           unsigned int j = context.dof_indices_var[_mesh_x_var][i];
00958           delta_x[i] = libmesh_real(this->current_solution(j)) -
00959                        libmesh_real(unsteady->old_nonlinear_solution(j));
00960         }
00961 
00962       for (unsigned int i = 0; i != n_y_dofs; ++i)
00963         {
00964           unsigned int j = context.dof_indices_var[_mesh_y_var][i];
00965           delta_y[i] = libmesh_real(this->current_solution(j)) -
00966                        libmesh_real(unsteady->old_nonlinear_solution(j));
00967         }
00968 
00969       for (unsigned int i = 0; i != n_z_dofs; ++i)
00970         {
00971           unsigned int j = context.dof_indices_var[_mesh_z_var][i];
00972           delta_z[i] = libmesh_real(this->current_solution(j)) -
00973                        libmesh_real(unsteady->old_nonlinear_solution(j));
00974         }
00975 
00976       for (unsigned int qp = 0; qp != n_qpoints; ++qp)
00977         {
00978           Gradient grad_u = context.interior_gradient(var, qp);
00979           RealGradient convection(0.);
00980 
00981           for (unsigned int i = 0; i != n_x_dofs; ++i)
00982             convection(0) += delta_x[i] * psi[i][qp];
00983           for (unsigned int i = 0; i != n_y_dofs; ++i)
00984             convection(1) += delta_y[i] * psi[i][qp];
00985           for (unsigned int i = 0; i != n_z_dofs; ++i)
00986             convection(2) += delta_z[i] * psi[i][qp];
00987 
00988           for (unsigned int i = 0; i != n_u_dofs; ++i)
00989             {
00990               Number JxWxPhiI = JxW[qp] * phi[i][qp];
00991               Fu(i) += (convection * grad_u) * JxWxPhiI;
00992               if (request_jacobian)
00993                 {
00994                   Number JxWxPhiI = JxW[qp] * phi[i][qp];
00995                   for (unsigned int j = 0; j != n_u_dofs; ++j)
00996                     Kuu(i,j) += JxWxPhiI * (convection * dphi[j][qp]);
00997 
00998                   Number JxWxPhiIoverDT = JxWxPhiI/this->deltat;
00999 
01000                   Number JxWxPhiIxDUDXoverDT = JxWxPhiIoverDT * grad_u(0);
01001                   for (unsigned int j = 0; j != n_x_dofs; ++j)
01002                     (*Kux)(i,j) += JxWxPhiIxDUDXoverDT * psi[j][qp];
01003 
01004                   Number JxWxPhiIxDUDYoverDT = JxWxPhiIoverDT * grad_u(1);
01005                   for (unsigned int j = 0; j != n_y_dofs; ++j)
01006                     (*Kuy)(i,j) += JxWxPhiIxDUDYoverDT * psi[j][qp];
01007 
01008                   Number JxWxPhiIxDUDZoverDT = JxWxPhiIoverDT * grad_u(2);
01009                   for (unsigned int j = 0; j != n_z_dofs; ++j)
01010                     (*Kuz)(i,j) += JxWxPhiIxDUDZoverDT * psi[j][qp];
01011                 }
01012             }
01013         }
01014     }
01015 
01016   return request_jacobian;
01017 }

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

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

Uses the forward sensitivity method.

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

Reimplemented from System.

Definition at line 749 of file implicit_system.C.

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

00752 {
00753   const unsigned int Np = parameters.size();
00754   const unsigned int Nq = qoi.size();
00755 
00756   // We currently get partial derivatives via central differencing
00757   const Real delta_p = TOLERANCE;
00758 
00759   // An introduction to the problem:
00760   //
00761   // Residual R(u(p),p) = 0
00762   // partial R / partial u = J = system matrix
00763   //
00764   // This implies that:
00765   // d/dp(R) = 0
00766   // (partial R / partial p) + 
00767   // (partial R / partial u) * (partial u / partial p) = 0
00768 
00769   // We first solve for (partial u / partial p) for each parameter:
00770   // J * (partial u / partial p) = - (partial R / partial p)
00771 
00772   this->sensitivity_solve(parameters);
00773 
00774   // Get ready to fill in senstivities:
00775   sensitivities.allocate_data(qoi_indices, *this, parameters);
00776 
00777   // We use the identity:
00778   // dq/dp = (partial q / partial p) + (partial q / partial u) *
00779   //         (partial u / partial p)
00780  
00781   // We get (partial q / partial u) from the user
00782   this->assemble_qoi_derivative(qoi_indices);
00783 
00784   for (unsigned int j=0; j != Np; ++j)
00785     {
00786       // (partial q / partial p) ~= (q(p+dp)-q(p-dp))/(2*dp)
00787 
00788       Number old_parameter = *parameters[j];
00789 
00790       *parameters[j] = old_parameter - delta_p;
00791       this->assemble_qoi();
00792       std::vector<Number> qoi_minus = this->qoi;
00793 
00794       *parameters[j] = old_parameter + delta_p;
00795       this->assemble_qoi();
00796       std::vector<Number>& qoi_plus = this->qoi;
00797 
00798       std::vector<Number> partialq_partialp(Nq, 0);
00799       for (unsigned int i=0; i != Nq; ++i)
00800         if (qoi_indices.has_index(i))
00801           partialq_partialp[i] = (qoi_plus[i] - qoi_minus[i]) / (2.*delta_p);
00802 
00803       // Don't leave the parameter changed
00804       *parameters[j] = old_parameter;
00805 
00806       for (unsigned int i=0; i != Nq; ++i)
00807         if (qoi_indices.has_index(i))
00808           sensitivities[i][j] = partialq_partialp[i] +
00809             this->get_adjoint_rhs(i).dot(this->get_sensitivity_solution(i));
00810     }
00811 
00812   // All parameters have been reset.
00813   // We didn't cache the original rhs or matrix for memory reasons,
00814   // but we can restore them to a state consistent solution -
00815   // principle of least surprise.
00816   this->assembly(true, true);
00817   this->rhs->close();
00818   this->matrix->close();
00819   this->assemble_qoi(qoi_indices);
00820 }

const NumericVector< Number > & System::get_adjoint_rhs ( unsigned int i = 0 ) const [inherited]

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

Definition at line 844 of file system.C.

References System::get_vector().

00845 {
00846   OStringStream adjoint_rhs_name;
00847   adjoint_rhs_name << "adjoint_rhs" << i;
00848 
00849   return this->get_vector(adjoint_rhs_name.str());
00850 }

NumericVector< Number > & System::get_adjoint_rhs ( unsigned int i = 0 ) [inherited]

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

Definition at line 834 of file system.C.

References System::get_vector().

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

00835 {
00836   OStringStream adjoint_rhs_name;
00837   adjoint_rhs_name << "adjoint_rhs" << i;
00838 
00839   return this->get_vector(adjoint_rhs_name.str());
00840 }

const NumericVector< Number > & System::get_adjoint_solution ( unsigned int i = 0 ) const [inherited]

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

Definition at line 784 of file system.C.

References System::get_vector().

00785 {
00786   OStringStream adjoint_name;
00787   adjoint_name << "adjoint_solution" << i;
00788 
00789   return this->get_vector(adjoint_name.str());
00790 }

NumericVector< Number > & System::get_adjoint_solution ( unsigned int i = 0 ) [inherited]

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

Definition at line 774 of file system.C.

References System::get_vector().

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

00775 {
00776   OStringStream adjoint_name;
00777   adjoint_name << "adjoint_solution" << i;
00778 
00779   return this->get_vector(adjoint_name.str());
00780 }

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

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

Definition at line 1364 of file system.h.

References System::_dof_map.

01365 {
01366   return *_dof_map;
01367 }

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

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

Definition at line 1356 of file system.h.

References System::_dof_map.

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

01357 {
01358   return *_dof_map;
01359 }

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

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

Definition at line 389 of file system.h.

References System::_equation_systems.

00389 { return _equation_systems; }

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

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

Definition at line 384 of file system.h.

References System::_equation_systems.

00384 { return _equation_systems; }

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

Gets a string containing the reference information.

Definition at line 45 of file reference_counter.C.

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

Referenced by ReferenceCounter::print_info().

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

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

Returns:: a string containing information about the system.

Definition at line 1233 of file system.C.

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

01234 {
01235   std::ostringstream out;
01236 
01237   
01238   const std::string& sys_name = this->name();
01239       
01240   out << "   System \"" << sys_name << "\"\n"
01241       << "    Type \""  << this->system_type() << "\"\n"
01242       << "    Variables=";
01243   
01244   for (unsigned int vn=0; vn<this->n_vars(); vn++)
01245       out << "\"" << this->variable_name(vn) << "\" ";
01246      
01247   out << '\n';
01248 
01249   out << "    Finite Element Types=";
01250 #ifndef LIBMESH_ENABLE_INFINITE_ELEMENTS
01251   for (unsigned int vn=0; vn<this->n_vars(); vn++)
01252     out << "\""
01253         << Utility::enum_to_string<FEFamily>(this->get_dof_map().variable_type(vn).family)
01254         << "\" ";  
01255 #else
01256   for (unsigned int vn=0; vn<this->n_vars(); vn++)
01257     {
01258       out << "\""
01259           << Utility::enum_to_string<FEFamily>(this->get_dof_map().variable_type(vn).family)
01260           << "\", \""
01261           << Utility::enum_to_string<FEFamily>(this->get_dof_map().variable_type(vn).radial_family)
01262           << "\" ";
01263     }
01264 
01265   out << '\n' << "    Infinite Element Mapping=";
01266   for (unsigned int vn=0; vn<this->n_vars(); vn++)
01267     out << "\""
01268         << Utility::enum_to_string<InfMapType>(this->get_dof_map().variable_type(vn).inf_map)
01269         << "\" ";
01270 #endif      
01271 
01272   out << '\n';
01273       
01274   out << "    Approximation Orders=";
01275   for (unsigned int vn=0; vn<this->n_vars(); vn++)
01276     {
01277 #ifndef LIBMESH_ENABLE_INFINITE_ELEMENTS
01278       out << "\""
01279           << Utility::enum_to_string<Order>(this->get_dof_map().variable_type(vn).order)
01280           << "\" ";
01281 #else
01282       out << "\""
01283           << Utility::enum_to_string<Order>(this->get_dof_map().variable_type(vn).order)
01284           << "\", \""
01285           << Utility::enum_to_string<Order>(this->get_dof_map().variable_type(vn).radial_order)
01286           << "\" ";
01287 #endif
01288     }
01289 
01290   out << '\n';
01291       
01292   out << "    n_dofs()="             << this->n_dofs()             << '\n';
01293   out << "    n_local_dofs()="       << this->n_local_dofs()       << '\n';
01294 #ifdef LIBMESH_ENABLE_AMR
01295   out << "    n_constrained_dofs()=" << this->n_constrained_dofs() << '\n';
01296 #endif
01297 
01298   out << "    " << "n_vectors()="  << this->n_vectors()  << '\n';
01299 //   out << "    " << "n_additional_matrices()=" << this->n_additional_matrices() << '\n';
01300   
01301   return out.str();
01302 }

std::pair< unsigned int, Real > DifferentiableSystem::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 from ImplicitSystem.

Definition at line 123 of file diff_system.C.

References DifferentiableSystem::time_solver.

00124 {
00125   return std::make_pair(this->time_solver->diff_solver()->max_linear_iterations,
00126                         this->time_solver->diff_solver()->relative_residual_tolerance);
00127 }

LinearSolver< Number > * DifferentiableSystem::get_linear_solver ( ) const [virtual, inherited]

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

Reimplemented from ImplicitSystem.

Definition at line 116 of file diff_system.C.

References AutoPtr< Tp >::get(), and DifferentiableSystem::time_solver.

00117 {
00118   return this->time_solver->linear_solver().get();
00119 }

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

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

Definition at line 245 of file implicit_system.C.

References ImplicitSystem::_matrices.

00246 {
00247   // Make sure the matrix exists
00248   matrices_iterator pos = _matrices.find (mat_name);
00249   
00250   if (pos == _matrices.end())
00251     {
00252       std::cerr << "ERROR: matrix "
00253                 << mat_name
00254                 << " does not exist in this system!"
00255                 << std::endl;      
00256       libmesh_error();
00257     }
00258   
00259   return *(pos->second);
00260 }

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

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

Definition at line 226 of file implicit_system.C.

References ImplicitSystem::_matrices.

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

00227 {
00228   // Make sure the matrix exists
00229   const_matrices_iterator pos = _matrices.find (mat_name);
00230   
00231   if (pos == _matrices.end())
00232     {
00233       std::cerr << "ERROR: matrix "
00234                 << mat_name
00235                 << " does not exist in this system!"
00236                 << std::endl;      
00237       libmesh_error();
00238     }
00239   
00240   return *(pos->second);
00241 }

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

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

Definition at line 1348 of file system.h.

References System::_mesh.

01349 {
01350   return _mesh;
01351 }

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

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

Definition at line 1340 of file system.h.

References System::_mesh.

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

01341 {
01342   return _mesh;
01343 }

const NumericVector< Number > & System::get_sensitivity_rhs ( unsigned int i = 0 ) const [inherited]

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

Definition at line 874 of file system.C.

References System::get_vector().

00875 {
00876   OStringStream sensitivity_rhs_name;
00877   sensitivity_rhs_name << "sensitivity_rhs" << i;
00878 
00879   return this->get_vector(sensitivity_rhs_name.str());
00880 }

NumericVector< Number > & System::get_sensitivity_rhs ( unsigned int i = 0 ) [inherited]

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

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

Definition at line 864 of file system.C.

References System::get_vector().

Referenced by ImplicitSystem::sensitivity_solve().

00865 {
00866   OStringStream sensitivity_rhs_name;
00867   sensitivity_rhs_name << "sensitivity_rhs" << i;
00868 
00869   return this->get_vector(sensitivity_rhs_name.str());
00870 }

const NumericVector< Number > & System::get_sensitivity_solution ( unsigned int i = 0 ) const [inherited]

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

Definition at line 733 of file system.C.

References System::get_vector().

00734 {
00735   OStringStream sensitivity_name;
00736   sensitivity_name << "sensitivity_solution" << i;
00737 
00738   return this->get_vector(sensitivity_name.str());
00739 }

NumericVector< Number > & System::get_sensitivity_solution ( unsigned int i = 0 ) [inherited]

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

Definition at line 723 of file system.C.

References System::get_vector().

Referenced by ImplicitSystem::sensitivity_solve().

00724 {
00725   OStringStream sensitivity_name;
00726   sensitivity_name << "sensitivity_solution" << i;
00727 
00728   return this->get_vector(sensitivity_name.str());
00729 }

NumericVector< Number > & System::get_vector ( const unsigned int vec_num ) [inherited]

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

Definition at line 682 of file system.C.

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

00683 {
00684   vectors_iterator v = vectors_begin();
00685   vectors_iterator v_end = vectors_end();
00686   unsigned int num = 0;
00687   while((num<vec_num) && (v!=v_end))
00688     {
00689       num++;
00690       ++v;
00691     }
00692   libmesh_assert(v!=v_end);
00693   return *(v->second);
00694 }

const NumericVector< Number > & System::get_vector ( const unsigned int vec_num ) const [inherited]

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

Definition at line 666 of file system.C.

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

00667 {
00668   const_vectors_iterator v = vectors_begin();
00669   const_vectors_iterator v_end = vectors_end();
00670   unsigned int num = 0;
00671   while((num<vec_num) && (v!=v_end))
00672     {
00673       num++;
00674       ++v;
00675     }
00676   libmesh_assert(v!=v_end);
00677   return *(v->second);
00678 }

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

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

Definition at line 647 of file system.C.

References System::_vectors.

00648 {
00649   // Make sure the vector exists
00650   vectors_iterator pos = _vectors.find(vec_name);
00651   
00652   if (pos == _vectors.end())
00653     {
00654       std::cerr << "ERROR: vector "
00655                 << vec_name
00656                 << " does not exist in this system!"
00657                 << std::endl;      
00658       libmesh_error();
00659     }
00660   
00661   return *(pos->second);
00662 }

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

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

Definition at line 628 of file system.C.

References System::_vectors.

Referenced by System::add_weighted_sensitivity_adjoint_solution(), UnsteadySolver::advance_timestep(), AdaptiveTimeSolver::advance_timestep(), System::compare(), UnsteadySolver::du(), System::get_adjoint_rhs(), System::get_adjoint_solution(), System::get_sensitivity_rhs(), System::get_sensitivity_solution(), System::get_weighted_sensitivity_adjoint_solution(), System::get_weighted_sensitivity_solution(), NewmarkSystem::initial_conditions(), UnsteadySolver::solve(), TwostepTimeSolver::solve(), FrequencySystem::solve(), NewmarkSystem::update_rhs(), and NewmarkSystem::update_u_v_a().

00629 {
00630   // Make sure the vector exists
00631   const_vectors_iterator pos = _vectors.find(vec_name);
00632   
00633   if (pos == _vectors.end())
00634     {
00635       std::cerr << "ERROR: vector "
00636                 << vec_name
00637                 << " does not exist in this system!"
00638                 << std::endl;      
00639       libmesh_error();
00640     }
00641   
00642   return *(pos->second);
00643 }

const NumericVector< Number > & System::get_weighted_sensitivity_adjoint_solution ( unsigned int i = 0 ) const [inherited]

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

Definition at line 814 of file system.C.

References System::get_vector().

00815 {
00816   OStringStream adjoint_name;
00817   adjoint_name << "weighted_sensitivity_adjoint_solution" << i;
00818 
00819   return this->get_vector(adjoint_name.str());
00820 }

NumericVector< Number > & System::get_weighted_sensitivity_adjoint_solution ( unsigned int i = 0 ) [inherited]

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

Definition at line 804 of file system.C.

References System::get_vector().

Referenced by ImplicitSystem::weighted_sensitivity_adjoint_solve().

00805 {
00806   OStringStream adjoint_name;
00807   adjoint_name << "weighted_sensitivity_adjoint_solution" << i;
00808 
00809   return this->get_vector(adjoint_name.str());
00810 }

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

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

Definition at line 757 of file system.C.

References System::get_vector().

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

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

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

Definition at line 750 of file system.C.

References System::get_vector().

Referenced by ImplicitSystem::weighted_sensitivity_solve().

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

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

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

Definition at line 935 of file system.C.

References System::_variable_numbers.

Referenced by GMVIO::copy_nodal_solution().

00936 {
00937   return _variable_numbers.count(var);
00938 }

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

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

Definition at line 369 of file implicit_system.h.

References ImplicitSystem::_matrices.

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

00370 {
00371   return (_matrices.count(mat_name));
00372 }

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

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

Definition at line 1450 of file system.h.

References System::_vectors.

Referenced by System::add_vector().

01451 {
01452   return (_vectors.count(vec_name));
01453 }

void ReferenceCounter::increment_constructor_count ( const std::string & name ) [inline, protected, inherited]

Increments the construction counter. Should be called in the constructor of any derived class that will be reference counted.

Definition at line 149 of file reference_counter.h.

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

Referenced by ReferenceCountedObject< SparseMatrix< T > >::ReferenceCountedObject().

00150 {
00151   Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
00152   std::pair<unsigned int, unsigned int>& p = _counts[name];
00153 
00154   p.first++;
00155 }

void ReferenceCounter::increment_destructor_count ( const std::string & name ) [inline, protected, inherited]

Increments the destruction counter. Should be called in the destructor of any derived class that will be reference counted.

Definition at line 167 of file reference_counter.h.

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

Referenced by ReferenceCountedObject< SparseMatrix< T > >::~ReferenceCountedObject().

00168 {
00169   Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
00170   std::pair<unsigned int, unsigned int>& p = _counts[name];
00171 
00172   p.second++;
00173 }

void System::init ( ) [inherited]

Initializes degrees of freedom on the current mesh. Sets the

Definition at line 158 of file system.C.

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

00159 { 
00160   // First initialize any required data
00161   this->init_data();
00162   
00163   //If no variables have been added to this system
00164   //don't do anything
00165   if(!this->n_vars())
00166     return;
00167  
00168   // Then call the user-provided intialization function
00169   this->user_initialization();
00170 }

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

Reimplemented from DifferentiableSystem.

Definition at line 749 of file fem_system.C.

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

Referenced by mesh_position_get().

00750 {
00751   Parent::init_context(c);
00752 
00753   FEMContext &context = libmesh_cast_ref<FEMContext&>(c);
00754 
00755   // Make sure we're prepared to do mass integration
00756   for (unsigned int var = 0; var != this->n_vars(); ++var)
00757     if (_time_evolving[var])
00758       {
00759         context.element_fe_var[var]->get_JxW();
00760         context.element_fe_var[var]->get_phi();
00761       }
00762 }

void FEMSystem::init_data ( ) [protected, virtual]

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

Reimplemented from DifferentiableSystem.

Reimplemented in ContinuationSystem.

Definition at line 71 of file fem_system.C.

References DifferentiableSystem::init_data().

Referenced by ContinuationSystem::init_data().

00072 {
00073   // First initialize LinearImplicitSystem data
00074   Parent::init_data();
00075 }

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

Initializes the matrices associated with this system.

Definition at line 99 of file implicit_system.C.

References ImplicitSystem::_can_add_matrices, ImplicitSystem::_matrices, DofMap::attach_matrix(), DofMap::compute_sparsity(), System::get_dof_map(), System::get_mesh(), and ImplicitSystem::matrix.

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

00100 {
00101   libmesh_assert (matrix != NULL);
00102 
00103   // Check for quick return in case the system matrix
00104   // (and by extension all the matrices) has already
00105   // been initialized
00106   if (matrix->initialized())
00107     return;
00108 
00109   // Get a reference to the DofMap
00110   DofMap& dof_map = this->get_dof_map();
00111   
00112   // no chance to add other matrices
00113   _can_add_matrices = false;
00114   
00115   // Tell the matrices about the dof map, and vice versa
00116   for (matrices_iterator pos = _matrices.begin();
00117        pos != _matrices.end(); ++pos)
00118     {
00119       libmesh_assert (!pos->second->initialized());
00120       dof_map.attach_matrix (*(pos->second));
00121     }
00122   
00123   // Compute the sparsity pattern for the current
00124   // mesh and DOF distribution.  This also updates
00125   // additional matrices, \p DofMap now knows them
00126   dof_map.compute_sparsity (this->get_mesh());
00127   
00128   // Initialize matrices
00129   for (matrices_iterator pos = _matrices.begin(); 
00130        pos != _matrices.end(); ++pos)
00131     pos->second->init ();
00132   
00133   // Set the additional matrices to 0.
00134   for (matrices_iterator pos = _matrices.begin(); 
00135        pos != _matrices.end(); ++pos)
00136     pos->second->zero ();
00137 }

void System::local_dof_indices	(	const unsigned int	var,
		std::set< unsigned int > &	var_indices
	)			const `[inherited]`

Fills the std::set with the degrees of freedom on the local processor corresponding the the variable number passed in.

Definition at line 962 of file system.C.

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

Referenced by System::discrete_var_norm().

00963 {
00964   // Make sure the set is clear
00965   var_indices.clear();
00966 
00967   std::vector<unsigned int> dof_indices;
00968   
00969   // Begin the loop over the elements
00970   MeshBase::const_element_iterator       el     =
00971     this->get_mesh().active_local_elements_begin();
00972   const MeshBase::const_element_iterator end_el =
00973     this->get_mesh().active_local_elements_end();
00974 
00975   const unsigned int 
00976     first_local = this->get_dof_map().first_dof(),
00977     end_local   = this->get_dof_map().end_dof();
00978 
00979   for ( ; el != end_el; ++el)
00980     {
00981       const Elem* elem = *el;      
00982       this->get_dof_map().dof_indices (elem, dof_indices, var);
00983 
00984       for(unsigned int i=0; i<dof_indices.size(); i++)
00985         {
00986           unsigned int dof = dof_indices[i];
00987 
00988           //If the dof is owned by the local processor
00989           if(first_local <= dof && dof < end_local)
00990             var_indices.insert(dof_indices[i]);
00991         }
00992     }
00993 }

bool FEMSystem::mass_residual	(	bool	request_jacobian,
		DiffContext &	context
	)			`[virtual]`

Adds a mass vector contribution on elem to elem_residual. If this method receives request_jacobian = true, then it should compute elem_jacobian and return true if possible. If elem_jacobian has not been computed then the method should return false.

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

Reimplemented from DifferentiableSystem.

Definition at line 1021 of file fem_system.C.

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

01023 {
01024   FEMContext &context = libmesh_cast_ref<FEMContext&>(c);
01025 
01026   unsigned int n_qpoints = context.element_qrule->n_points();
01027 
01028   for (unsigned int var = 0; var != this->n_vars(); ++var)
01029     {
01030       if (!_time_evolving[var])
01031         continue;
01032 
01033       const std::vector<Real> &JxW = 
01034         context.element_fe_var[var]->get_JxW();
01035 
01036       const std::vector<std::vector<Real> > &phi =
01037         context.element_fe_var[var]->get_phi();
01038 
01039       const unsigned int n_dofs = context.dof_indices_var[var].size();
01040 
01041       DenseSubVector<Number> &Fu = *context.elem_subresiduals[var];
01042       DenseSubMatrix<Number> &Kuu = *context.elem_subjacobians[var][var];
01043 
01044       for (unsigned int qp = 0; qp != n_qpoints; ++qp)
01045         {
01046           Number u = context.interior_value(var, qp);
01047           Number JxWxU = JxW[qp] * u;
01048           for (unsigned int i = 0; i != n_dofs; ++i)
01049             {
01050               Fu(i) += JxWxU * phi[i][qp];
01051               if (request_jacobian && context.elem_solution_derivative)
01052                 {
01053                   libmesh_assert (context.elem_solution_derivative == 1.0);
01054 
01055                   Number JxWxPhiI = JxW[qp] * phi[i][qp];
01056                   Kuu(i,i) += JxWxPhiI * phi[i][qp];
01057                   for (unsigned int j = i+1; j < n_dofs; ++j)
01058                     {
01059                       Number Kij = JxWxPhiI * phi[j][qp];
01060                       Kuu(i,j) += Kij;
01061                       Kuu(j,i) += Kij;
01062                     }
01063                 }
01064             }
01065         }
01066     }
01067 
01068   return request_jacobian;
01069 }

void FEMSystem::mesh_position_get ( )

Tells the FEMSystem to set the degree of freedom coefficients which should correspond to mesh nodal coordinates.

Definition at line 810 of file fem_system.C.

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

00811 {
00812   // This function makes no sense unless we've already picked out some
00813   // variable(s) to reflect mesh position coordinates
00814   if (_mesh_sys == libMesh::invalid_uint)
00815     libmesh_error();
00816 
00817   // We currently assume mesh variables are in our own system
00818   if (_mesh_sys != this->number())
00819     libmesh_not_implemented();
00820 
00821   // Loop over every active mesh element on this processor
00822   const MeshBase& mesh = this->get_mesh();
00823 
00824   const DofMap& dof_map = this->get_dof_map();
00825 
00826   MeshBase::const_element_iterator el =
00827     mesh.active_local_elements_begin();
00828   const MeshBase::const_element_iterator end_el =
00829     mesh.active_local_elements_end();
00830 
00831   AutoPtr<DiffContext> con = this->build_context();
00832   FEMContext &_femcontext = libmesh_cast_ref<FEMContext&>(*con);
00833   this->init_context(_femcontext);
00834 
00835   // Get the solution's mesh variables from every element
00836   for ( ; el != end_el; ++el)
00837     {
00838       _femcontext.elem = *el;
00839 
00840       // Initialize the per-variable data for elem.
00841       for (unsigned int i=0; i != this->n_vars(); ++i)
00842         {
00843           dof_map.dof_indices (_femcontext.elem,
00844                                _femcontext.dof_indices_var[i], i);
00845         }
00846 
00847       _femcontext.elem_position_get();
00848       if (_mesh_x_var != libMesh::invalid_uint)
00849         this->solution->insert(*_femcontext.elem_subsolutions[_mesh_x_var],
00850                                _femcontext.dof_indices_var[_mesh_x_var]);
00851       if (_mesh_y_var != libMesh::invalid_uint)
00852         this->solution->insert(*_femcontext.elem_subsolutions[_mesh_y_var],
00853                                _femcontext.dof_indices_var[_mesh_y_var]);
00854       if (_mesh_z_var != libMesh::invalid_uint)
00855         this->solution->insert(*_femcontext.elem_subsolutions[_mesh_z_var],
00856                                _femcontext.dof_indices_var[_mesh_z_var]);
00857     }
00858 
00859   // And make sure the current_local_solution is up to date too
00860   this->update();
00861 }

void FEMSystem::mesh_position_set ( )

Tells the FEMSystem to set the mesh nodal coordinates which should correspond to degree of freedom coefficients.

Definition at line 392 of file fem_system.C.

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

Referenced by solve().

00393 {
00394   // If we don't need to move the mesh, we're done
00395   if (_mesh_sys != this->number())
00396     return;
00397 
00398   const MeshBase& mesh = this->get_mesh();
00399 
00400   // This code won't work on a parallelized mesh yet -
00401   // it won't get ancestor elements right.
00402   libmesh_assert(mesh.is_serial());
00403   
00404   AutoPtr<DiffContext> con = this->build_context();
00405   FEMContext &_femcontext = libmesh_cast_ref<FEMContext&>(*con);
00406 
00407   // Move every mesh element we can
00408   MeshBase::const_element_iterator el =
00409     mesh.active_elements_begin();
00410   const MeshBase::const_element_iterator end_el =
00411     mesh.active_elements_end();
00412 
00413   for ( ; el != end_el; ++el)
00414     {
00415       _femcontext.reinit(*this, *el);
00416 
00417       // This code won't handle moving subactive elements
00418       libmesh_assert(!_femcontext.elem->has_children());
00419 
00420       _femcontext.elem_position_set(0.);
00421     }
00422 }

void FEMSystem::mesh_x_position	(	unsigned int	sysnum,
		unsigned int	var
	)			`[virtual]`

Tells the FEMSystem that variable var from system number sysnum should be used to update the x coordinate of mesh nodes, in problems where the mesh itself is expected to move in time.

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

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

Definition at line 774 of file fem_system.C.

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

00775 {
00776   if (_mesh_sys != libMesh::invalid_uint && _mesh_sys != sysnum)
00777     libmesh_error();
00778   if (sysnum != this->number())
00779     libmesh_not_implemented();
00780   _mesh_sys = sysnum;
00781   _mesh_x_var = var;
00782 }

void FEMSystem::mesh_y_position	(	unsigned int	sysnum,
		unsigned int	var
	)			`[virtual]`

Tells the FEMSystem that variable var from system number sysnum should be used to update the y coordinate of mesh nodes.

Definition at line 786 of file fem_system.C.

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

00787 {
00788   if (_mesh_sys != libMesh::invalid_uint && _mesh_sys != sysnum)
00789     libmesh_error();
00790   if (sysnum != this->number())
00791     libmesh_not_implemented();
00792   _mesh_sys = sysnum;
00793   _mesh_y_var = var;
00794 }

void FEMSystem::mesh_z_position	(	unsigned int	sysnum,
		unsigned int	var
	)			`[virtual]`

Tells the FEMSystem that variable var from system number sysnum should be used to update the z coordinate of mesh nodes.

Definition at line 798 of file fem_system.C.

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

00799 {
00800   if (_mesh_sys != libMesh::invalid_uint && _mesh_sys != sysnum)
00801     libmesh_error();
00802   if (sysnum != this->number())
00803     libmesh_not_implemented();
00804   _mesh_sys = sysnum;
00805   _mesh_z_var = var;
00806 }

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

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

Definition at line 1442 of file system.h.

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

01443 {
01444   return this->n_dofs() - this->n_constrained_dofs();
01445 }

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

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

Definition at line 95 of file system.C.

References System::_dof_map.

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

00096 {
00097 #ifdef LIBMESH_ENABLE_AMR
00098 
00099   return _dof_map->n_constrained_dofs();
00100 
00101 #else
00102 
00103   return 0;
00104 
00105 #endif
00106 }

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

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

Definition at line 88 of file system.C.

References System::_dof_map.

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

00089 {
00090   return _dof_map->n_dofs();
00091 }

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

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

Definition at line 110 of file system.C.

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

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

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

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

Returns:: the number of matrices handled by this system

Definition at line 376 of file implicit_system.h.

References ImplicitSystem::_matrices.

00377 {
00378  return _matrices.size(); 
00379 }

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

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

Definition at line 76 of file reference_counter.h.

References ReferenceCounter::_n_objects.

00077   { return _n_objects; }

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

Returns:: the number of variables in the system

Definition at line 1396 of file system.h.

References System::_variables.

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

01397 {
01398   return _variables.size();
01399 }

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

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

Definition at line 1458 of file system.h.

References System::_vectors.

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

01459 {
01460   return _vectors.size(); 
01461 }

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

Returns:: the system name.

Definition at line 1324 of file system.h.

References System::_sys_name.

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

01325 {
01326   return _sys_name;
01327 }

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

Returns:: the system number.

Definition at line 1332 of file system.h.

References System::_sys_number.

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

01333 {
01334   return _sys_number;
01335 }

void FEMSystem::numerical_elem_jacobian ( FEMContext & context ) [protected]

Uses the results of multiple element_residual() calls to numerically differentiate the corresponding jacobian on an element.

Definition at line 724 of file fem_system.C.

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

Referenced by assembly().

00725 {
00726   START_LOG("numerical_elem_jacobian()", "FEMSystem");
00727   this->numerical_jacobian(&TimeSolver::element_residual, context);
00728   STOP_LOG("numerical_elem_jacobian()", "FEMSystem");
00729 }

void FEMSystem::numerical_jacobian	(	TimeSolverResPtr	res,
		FEMContext &	context
	)			`[protected]`

Uses the results of multiple res calls to numerically differentiate the corresponding jacobian.

Definition at line 621 of file fem_system.C.

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

Referenced by numerical_elem_jacobian(), and numerical_side_jacobian().

00623 {
00624   // Logging is done by numerical_elem_jacobian
00625   // or numerical_side_jacobian
00626 
00627   DenseVector<Number> original_residual(context.elem_residual);
00628   DenseVector<Number> backwards_residual(context.elem_residual);
00629   DenseMatrix<Number> numerical_jacobian(context.elem_jacobian);
00630 #ifdef DEBUG
00631   DenseMatrix<Number> old_jacobian(context.elem_jacobian);
00632 #endif
00633 
00634   Real numerical_point_h = 0.;
00635   if (_mesh_sys == this->number())
00636     numerical_point_h = numerical_jacobian_h * context.elem->hmin();
00637 
00638   for (unsigned int j = 0; j != context.dof_indices.size(); ++j)
00639     {
00640       // Take the "minus" side of a central differenced first derivative
00641       Number original_solution = context.elem_solution(j);
00642       context.elem_solution(j) -= numerical_jacobian_h;
00643 
00644       // Make sure to catch any moving mesh terms
00645       // FIXME - this could be less ugly
00646       Real *coord = NULL;
00647       if (_mesh_sys == this->number())
00648         {
00649           if (_mesh_x_var != libMesh::invalid_uint)
00650             for (unsigned int k = 0;
00651                  k != context.dof_indices_var[_mesh_x_var].size(); ++k)
00652               if (context.dof_indices_var[_mesh_x_var][k] ==
00653                   context.dof_indices[j])
00654                 coord = &(context.elem->point(k)(0));
00655           if (_mesh_y_var != libMesh::invalid_uint)
00656             for (unsigned int k = 0;
00657                  k != context.dof_indices_var[_mesh_y_var].size(); ++k)
00658               if (context.dof_indices_var[_mesh_y_var][k] == 
00659                   context.dof_indices[j])
00660                 coord = &(context.elem->point(k)(1));
00661           if (_mesh_z_var != libMesh::invalid_uint)
00662             for (unsigned int k = 0;
00663                  k != context.dof_indices_var[_mesh_z_var].size(); ++k)
00664               if (context.dof_indices_var[_mesh_z_var][k] ==
00665                   context.dof_indices[j])
00666                 coord = &(context.elem->point(k)(2));
00667         }
00668       if (coord)
00669         {
00670           // We have enough information to scale the perturbations
00671           // here appropriately
00672           context.elem_solution(j) = original_solution - numerical_point_h;
00673           *coord = libmesh_real(context.elem_solution(j));
00674         }
00675 
00676       context.elem_residual.zero();
00677       ((*time_solver).*(res))(false, context);
00678 #ifdef DEBUG
00679       libmesh_assert(old_jacobian == context.elem_jacobian);
00680 #endif
00681       backwards_residual = context.elem_residual;
00682 
00683       // Take the "plus" side of a central differenced first derivative
00684       context.elem_solution(j) = original_solution + numerical_jacobian_h;
00685       if (coord)
00686         {
00687           context.elem_solution(j) = original_solution + numerical_point_h;
00688           *coord = libmesh_real(context.elem_solution(j));
00689         }
00690       context.elem_residual.zero();
00691       ((*time_solver).*(res))(false, context);
00692 #ifdef DEBUG
00693       libmesh_assert(old_jacobian == context.elem_jacobian);
00694 #endif
00695 
00696       context.elem_solution(j) = original_solution;
00697       if (coord)
00698         {
00699           *coord = libmesh_real(context.elem_solution(j));
00700           for (unsigned int i = 0; i != context.dof_indices.size(); ++i)
00701             {
00702               numerical_jacobian(i,j) =
00703                 (context.elem_residual(i) - backwards_residual(i)) /
00704                 2. / numerical_point_h;
00705             }
00706         }
00707       else
00708         {
00709           for (unsigned int i = 0; i != context.dof_indices.size(); ++i)
00710             {
00711               numerical_jacobian(i,j) =
00712                 (context.elem_residual(i) - backwards_residual(i)) /
00713                 2. / numerical_jacobian_h;
00714             }
00715         }
00716     }
00717 
00718   context.elem_residual = original_residual;
00719   context.elem_jacobian = numerical_jacobian;
00720 }

void FEMSystem::numerical_side_jacobian ( FEMContext & context ) [protected]

Uses the results of multiple side_residual() calls to numerically differentiate the corresponding jacobian on an element's side.

Definition at line 733 of file fem_system.C.

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

Referenced by assembly().

00734 {
00735   START_LOG("numerical_side_jacobian()", "FEMSystem");
00736   this->numerical_jacobian(&TimeSolver::side_residual, context);
00737   STOP_LOG("numerical_side_jacobian()", "FEMSystem");
00738 }

void FEMSystem::postprocess ( ) [virtual]

Runs a postprocessing loop over all elements, and if postprocess_sides is true over all sides.

Implements DifferentiableSystem.

Definition at line 426 of file fem_system.C.

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

00427 {
00428   START_LOG("postprocess()", "FEMSystem");
00429 
00430   const MeshBase& mesh = this->get_mesh();
00431 
00432   this->update();
00433 
00434   AutoPtr<DiffContext> con = this->build_context();
00435   FEMContext &_femcontext = libmesh_cast_ref<FEMContext&>(*con);
00436 
00437   // Loop over every active mesh element on this processor
00438   MeshBase::const_element_iterator el =
00439     mesh.active_local_elements_begin();
00440   const MeshBase::const_element_iterator end_el =
00441     mesh.active_local_elements_end();
00442 
00443   for ( ; el != end_el; ++el)
00444     {
00445       _femcontext.reinit(*this, *el);
00446 
00447       // Optionally initialize all the interior FE objects on elem.
00448       // if (fe_reinit_during_postprocess)
00449       // _femcontext.elem_fe_reinit();
00450       
00451       this->element_postprocess(_femcontext);
00452 
00453       for (_femcontext.side = 0;
00454            _femcontext.side != _femcontext.elem->n_sides();
00455            ++_femcontext.side)
00456         {
00457           // Don't compute on non-boundary sides unless requested
00458           if (!postprocess_sides ||
00459               (!compute_internal_sides &&
00460                _femcontext.elem->neighbor(_femcontext.side) != NULL))
00461             continue;
00462 
00463           // Optionally initialize all the interior FE objects on elem/side.
00464           // Logging of FE::reinit is done in the FE functions
00465           if (fe_reinit_during_postprocess)
00466             {
00467               std::map<FEType, FEBase *>::iterator fe_end =
00468                 _femcontext.side_fe.end();
00469               for (std::map<FEType, FEBase *>::iterator i =
00470                 _femcontext.side_fe.begin();
00471                    i != fe_end; ++i)
00472                 {
00473                   i->second->reinit(_femcontext.elem, _femcontext.side);
00474                 }
00475             }
00476 
00477           this->side_postprocess(_femcontext);
00478         }
00479     }
00480 
00481   STOP_LOG("postprocess()", "FEMSystem");
00482 }

void ReferenceCounter::print_info ( ) [static, inherited]

Prints the reference information to std::cout.

Definition at line 83 of file reference_counter.C.

References ReferenceCounter::get_info().

00084 {
00085 #if defined(LIBMESH_ENABLE_REFERENCE_COUNTING) && defined(DEBUG)
00086   
00087   std::cout << ReferenceCounter::get_info();
00088   
00089 #endif
00090 }

void System::project_solution	(	Number	fptrconst Point &p,const Parameters &parameters,const std::string &sys_name,const std::string &unknown_name,
		Gradient	gptrconst Point &p,const Parameters &parameters,const std::string &sys_name,const std::string &unknown_name,
		Parameters &	parameters
	)			const `[inherited]`

Projects the continuous functions onto the current solution.

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

Definition at line 237 of file system_projection.C.

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

00246 {
00247   this->project_vector(fptr, gptr, parameters, *solution);
00248 
00249   solution->localize(*current_local_solution);
00250 }

bool& System::project_solution_on_reinit ( void ) [inline, inherited]

Tells the System whether or not to project the solution vector onto new grids when the system is reinitialized. The solution will be projected unless project_solution_on_reinit() = false is called.

Definition at line 451 of file system.h.

References System::_solution_projection.

00452     { return _solution_projection; }

void System::project_vector	(	const NumericVector< Number > &	old_v,
		NumericVector< Number > &	new_v
	)			const `[protected, inherited]`

Projects the vector defined on the old mesh onto the new mesh. The original vector is unchanged and the new vector is passed through the second argument.

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

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

Definition at line 60 of file system_projection.C.

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

00062 {
00063   START_LOG ("project_vector()", "System");
00064 
00071   new_v.clear();
00072 
00073 #ifdef LIBMESH_ENABLE_AMR
00074 
00075   // Resize the new vector and get a serial version.
00076   NumericVector<Number> *new_vector_ptr = NULL;
00077   AutoPtr<NumericVector<Number> > new_vector_built;
00078   NumericVector<Number> *local_old_vector;
00079   AutoPtr<NumericVector<Number> > local_old_vector_built;
00080   const NumericVector<Number> *old_vector_ptr = NULL;
00081 
00082   ConstElemRange active_local_elem_range 
00083     (this->get_mesh().active_local_elements_begin(),
00084      this->get_mesh().active_local_elements_end());
00085 
00086   // If the old vector was uniprocessor, make the new
00087   // vector uniprocessor
00088   if (old_v.type() == SERIAL)
00089     {
00090       new_v.init (this->n_dofs(), false, SERIAL);
00091       new_vector_ptr = &new_v;
00092       old_vector_ptr = &old_v;
00093     }
00094 
00095   // Otherwise it is a parallel, distributed vector, which
00096   // we need to localize.
00097   else if (old_v.type() == PARALLEL)
00098     {
00099       // Build a send list for efficient localization
00100       BuildProjectionList projection_list(*this);
00101       Threads::parallel_reduce (active_local_elem_range, 
00102                                 projection_list);
00103 
00104       // Create a sorted, unique send_list
00105       projection_list.unique();
00106 
00107       new_v.init (this->n_dofs(), this->n_local_dofs(), false, PARALLEL);
00108       new_vector_built = NumericVector<Number>::build();
00109       local_old_vector_built = NumericVector<Number>::build();
00110       new_vector_ptr = new_vector_built.get();
00111       local_old_vector = local_old_vector_built.get();
00112       new_vector_ptr->init(this->n_dofs(), false, SERIAL);
00113       local_old_vector->init(old_v.size(), false, SERIAL);
00114       old_v.localize(*local_old_vector, projection_list.send_list);
00115       local_old_vector->close();
00116       old_vector_ptr = local_old_vector;
00117     }
00118   else if (old_v.type() == GHOSTED)
00119     {
00120       // Build a send list for efficient localization
00121       BuildProjectionList projection_list(*this);
00122       Threads::parallel_reduce (active_local_elem_range, 
00123                                 projection_list);
00124 
00125       // Create a sorted, unique send_list
00126       projection_list.unique();
00127 
00128       new_v.init (this->n_dofs(), this->n_local_dofs(),
00129                   this->get_dof_map().get_send_list(), false, GHOSTED);
00130 
00131       local_old_vector_built = NumericVector<Number>::build();
00132       new_vector_ptr = &new_v;
00133       local_old_vector = local_old_vector_built.get();
00134       local_old_vector->init(old_v.size(), old_v.local_size(),
00135                              projection_list.send_list, false, GHOSTED);
00136       old_v.localize(*local_old_vector, projection_list.send_list);
00137       local_old_vector->close();
00138       old_vector_ptr = local_old_vector;
00139     }
00140   else // unknown old_v.type()
00141     {
00142       std::cerr << "ERROR: Unknown old_v.type() == " << old_v.type() 
00143                 << std::endl;
00144       libmesh_error();
00145     }
00146   
00147   // Note that the above will have zeroed the new_vector.
00148   // Just to be sure, assert that new_vector_ptr and old_vector_ptr
00149   // were successfully set before trying to deref them.
00150   libmesh_assert(new_vector_ptr);
00151   libmesh_assert(old_vector_ptr);
00152 
00153   NumericVector<Number> &new_vector = *new_vector_ptr;
00154   const NumericVector<Number> &old_vector = *old_vector_ptr;
00155     
00156   Threads::parallel_for (active_local_elem_range,
00157                          ProjectVector(*this,
00158                                        old_vector,
00159                                        new_vector)
00160                          );
00161 
00162   // Copy the SCALAR dofs from old_vector to new_vector
00163   // Note: We assume that all SCALAR dofs are on the
00164   // processor with highest ID
00165   if(libMesh::processor_id() == (libMesh::n_processors()-1))
00166   {
00167     const DofMap& dof_map = this->get_dof_map();
00168     for (unsigned int var=0; var<this->n_vars(); var++)
00169       if(this->variable(var).type().family == SCALAR)
00170         {
00171           // We can just map SCALAR dofs directly across
00172           std::vector<unsigned int> new_SCALAR_indices, old_SCALAR_indices;
00173           dof_map.SCALAR_dof_indices (new_SCALAR_indices, var, false);
00174           dof_map.SCALAR_dof_indices (old_SCALAR_indices, var, true);
00175           const unsigned int new_n_dofs = new_SCALAR_indices.size();
00176 
00177           for (unsigned int i=0; i<new_n_dofs; i++)
00178           {
00179             new_vector.set( new_SCALAR_indices[i], old_vector(old_SCALAR_indices[i]) );
00180           }
00181         }
00182   }
00183 
00184   new_vector.close();
00185 
00186   // If the old vector was serial, we probably need to send our values
00187   // to other processors
00188   //
00189   // FIXME: I'm not sure how to make a NumericVector do that without
00190   // creating a temporary parallel vector to use localize! - RHS
00191   if (old_v.type() == SERIAL)
00192     {
00193       AutoPtr<NumericVector<Number> > dist_v = NumericVector<Number>::build();
00194       dist_v->init(this->n_dofs(), this->n_local_dofs(), false, PARALLEL);
00195       dist_v->close();
00196     
00197       for (unsigned int i=0; i!=dist_v->size(); i++)
00198         if (new_vector(i) != 0.0)
00199           dist_v->set(i, new_vector(i));
00200 
00201       dist_v->close();
00202 
00203       dist_v->localize (new_v, this->get_dof_map().get_send_list());
00204       new_v.close();
00205     }
00206   // If the old vector was parallel, we need to update it
00207   // and free the localized copies
00208   else if (old_v.type() == PARALLEL)
00209     {
00210       // We may have to set dof values that this processor doesn't
00211       // own in certain special cases, like LAGRANGE FIRST or
00212       // HERMITE THIRD elements on second-order meshes
00213       for (unsigned int i=0; i!=new_v.size(); i++)
00214         if (new_vector(i) != 0.0)
00215           new_v.set(i, new_vector(i));
00216       new_v.close();
00217     }
00218 
00219   this->get_dof_map().enforce_constraints_exactly(*this, &new_v);
00220 
00221 #else
00222 
00223   // AMR is disabled: simply copy the vector
00224   new_v = old_v;
00225   
00226 #endif // #ifdef LIBMESH_ENABLE_AMR
00227 
00228   STOP_LOG("project_vector()", "System");
00229 }

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

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

Definition at line 43 of file system_projection.C.

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

00044 {
00045   // Create a copy of the vector, which currently
00046   // contains the old data.
00047   AutoPtr<NumericVector<Number> >
00048     old_vector (vector.clone());
00049 
00050   // Project the old vector to the new vector
00051   this->project_vector (*old_vector, vector);
00052 }

void System::project_vector	(	Number	fptrconst Point &p,const Parameters &parameters,const std::string &sys_name,const std::string &unknown_name,
		Gradient	gptrconst Point &p,const Parameters &parameters,const std::string &sys_name,const std::string &unknown_name,
		Parameters &	parameters,
		NumericVector< Number > &	new_vector
	)			const `[inherited]`

Projects the continuous functions onto the current mesh.

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

Definition at line 258 of file system_projection.C.

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

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

00268 {
00269   START_LOG ("project_vector()", "System");
00270 
00271   Threads::parallel_for (ConstElemRange (this->get_mesh().active_local_elements_begin(),
00272                                          this->get_mesh().active_local_elements_end(),
00273                                          1000),
00274                          ProjectSolution(*this,
00275                                          fptr,
00276                                          gptr,
00277                                          parameters,
00278                                          new_vector)
00279                          );
00280 
00281   // Also, load values into the SCALAR dofs
00282   // Note: We assume that all SCALAR dofs are on the
00283   // processor with highest ID
00284   if(libMesh::processor_id() == (libMesh::n_processors()-1))
00285   {
00286     const DofMap& dof_map = this->get_dof_map();
00287     for (unsigned int var=0; var<this->n_vars(); var++)
00288       if(this->variable(var).type().family == SCALAR)
00289         {
00290           std::vector<unsigned int> SCALAR_indices;
00291           dof_map.SCALAR_dof_indices (SCALAR_indices, var);
00292           const unsigned int n_SCALAR_dofs = SCALAR_indices.size();
00293 
00294           for (unsigned int i=0; i<n_SCALAR_dofs; i++)
00295           {
00296             const unsigned int index = SCALAR_indices[i];
00297 
00298             // We pass the point (i,0,0) to the fptr to distinguish
00299             // the different scalars within the SCALAR variable
00300             Point p_i(i,0,0);
00301             new_vector.set( index, fptr(p_i,
00302                                         parameters,
00303                                         this->name(),
00304                                         this->variable_name(var))
00305                           );
00306           }
00307         }
00308   }
00309 
00310   new_vector.close();
00311 
00312 #ifdef LIBMESH_ENABLE_AMR
00313   this->get_dof_map().enforce_constraints_exactly(*this, &new_vector);
00314 #endif
00315 
00316   STOP_LOG("project_vector()", "System");
00317 }

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

Prolong vectors after the mesh has refined

Definition at line 255 of file system.C.

References System::restrict_vectors().

Referenced by EquationSystems::reinit().

00256 {
00257 #ifdef LIBMESH_ENABLE_AMR
00258   // Currently project_vector handles both restriction and prolongation
00259   this->restrict_vectors();
00260 #endif
00261 }

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

Definition at line 1004 of file implicit_system.C.

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

01007 {
01008   // We currently get partial derivatives via finite differencing
01009   const Real delta_p = TOLERANCE;
01010 
01011   // We'll use one temporary vector for matrix-vector-vector products
01012   AutoPtr<NumericVector<Number> > tempvec = this->solution->zero_clone();
01013 
01014   // And another temporary vector to hold a copy of the true solution
01015   // so we can safely perturb this->solution.
01016   AutoPtr<NumericVector<Number> > oldsolution = this->solution->clone();
01017 
01018   const unsigned int Np = parameters.size();
01019   const unsigned int Nq = qoi.size();
01020 
01021   // For each quantity of interest q, the parameter sensitivity
01022   // Hessian is defined as q''_{kl} = {d^2 q}/{d p_k d p_l}.
01023   //
01024   // We calculate it from values and partial derivatives of the
01025   // quantity of interest function Q, solution u, adjoint solution z,
01026   // and residual R, as:
01027   //
01028   // q''_{kl} =
01029   // Q''_{kl} + Q''_{uk}(u)*u'_l + Q''_{ul}(u) * u'_k + 
01030   // Q''_{uu}(u)*u'_k*u'_l -
01031   // R''_{kl}(u,z) -
01032   // R''_{uk}(u,z)*u'_l - R''_{ul}(u,z)*u'_k -
01033   // R''_{uu}(u,z)*u'_k*u'_l
01034   //
01035   // See the adjoints model document for more details.
01036 
01037   // We first do an adjoint solve to get z for each quantity of
01038   // interest
01039 
01040   this->adjoint_solve(qoi_indices);
01041 
01042   // And a sensitivity solve to get u_k for each parameter
01043   this->sensitivity_solve(parameters);
01044 
01045   // Get ready to fill in second derivatives:
01046   sensitivities.allocate_hessian_data(qoi_indices, *this, parameters);
01047 
01048   for (unsigned int k=0; k != Np; ++k)
01049     {
01050       Number old_parameterk = *parameters[k];
01051 
01052       // The Hessian is symmetric, so we just calculate the lower
01053       // triangle and the diagonal, and we get the upper triangle from
01054       // the transpose of the lower
01055 
01056       for (unsigned int l=0; l != k+1; ++l)
01057         {
01058           // The second partial derivatives with respect to parameters
01059           // are all calculated via a central finite difference
01060           // stencil:
01061           // F''_{kl} ~= (F(p+dp*e_k+dp*e_l) - F(p+dp*e_k-dp*e_l) -
01062           //              F(p-dp*e_k+dp*e_l) - F(p-dp*e_k-dp*e_l))/(4*dp^2)
01063           // We will add Q''_{kl}(u) and subtract R''_{kl}(u,z) at the
01064           // same time.
01065 
01066           Number old_parameterl = *parameters[l];
01067 
01068           *parameters[k] = old_parameterk + delta_p;
01069           *parameters[l] = old_parameterl + delta_p;
01070           this->assemble_qoi(qoi_indices);
01071           this->assembly(true, false);
01072           std::vector<Number> partial2q_term = this->qoi;
01073           std::vector<Number> partial2R_term(this->qoi.size());
01074           for (unsigned int i=0; i != Nq; ++i)
01075             if (qoi_indices.has_index(i))
01076               partial2R_term[i] = this->rhs->dot(this->get_adjoint_solution(i));
01077 
01078           *parameters[l] = old_parameterl - delta_p;
01079           this->assemble_qoi(qoi_indices);
01080           this->assembly(true, false);
01081           for (unsigned int i=0; i != Nq; ++i)
01082             if (qoi_indices.has_index(i))
01083               {
01084                 partial2q_term[i] -= this->qoi[i];
01085                 partial2R_term[i] -= this->rhs->dot(this->get_adjoint_solution(i));
01086               }
01087 
01088           *parameters[k] = old_parameterk - delta_p;
01089           this->assemble_qoi(qoi_indices);
01090           this->assembly(true, false);
01091           for (unsigned int i=0; i != Nq; ++i)
01092             if (qoi_indices.has_index(i))
01093               {
01094                 partial2q_term[i] += this->qoi[i];
01095                 partial2R_term[i] += this->rhs->dot(this->get_adjoint_solution(i));
01096               }
01097 
01098           *parameters[l] = old_parameterl + delta_p;
01099           this->assemble_qoi(qoi_indices);
01100           this->assembly(true, false);
01101           for (unsigned int i=0; i != Nq; ++i)
01102             if (qoi_indices.has_index(i))
01103               {
01104                 partial2q_term[i] -= this->qoi[i];
01105                 partial2R_term[i] -= this->rhs->dot(this->get_adjoint_solution(i));
01106                 partial2q_term[i] /= (4. * delta_p * delta_p);
01107                 partial2R_term[i] /= (4. * delta_p * delta_p);
01108               }
01109 
01110           for (unsigned int i=0; i != Nq; ++i)
01111             if (qoi_indices.has_index(i))
01112               {
01113                 Number current_terms = partial2q_term[i] - partial2R_term[i];
01114                 sensitivities.second_derivative(i,k,l) += current_terms;
01115                 if (k != l)
01116                   sensitivities.second_derivative(i,l,k) += current_terms;
01117               }
01118 
01119           // Don't leave the parameters perturbed
01120           *parameters[l] = old_parameterl;
01121           *parameters[k] = old_parameterk;
01122         }
01123 
01124       // We get (partial q / partial u) and 
01125       // (partial R / partial u) from the user, but centrally
01126       // difference to get q_uk and R_uk terms:
01127       // (partial^2 q / partial u partial k)
01128       // q_uk*u'_l = (q_u(p+dp*e_k)*u'_l - q_u(p-dp*e_k)*u'_l)/(2*dp)
01129       // 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)
01130       //
01131       // To avoid creating Nq temporary vectors, we add these
01132       // subterms to the sensitivities output one by one.
01133       //
01134       // FIXME: this is probably a bad order of operations for
01135       // controlling floating point error.
01136 
01137       *parameters[k] = old_parameterk + delta_p;
01138       this->assembly(false, true);
01139       this->assemble_qoi_derivative(qoi_indices);
01140 
01141       for (unsigned int l=0; l != Np; ++l)
01142         {
01143           this->matrix->vector_mult(*tempvec, this->get_sensitivity_solution(l));
01144           for (unsigned int i=0; i != Nq; ++i)
01145             if (qoi_indices.has_index(i))
01146               {
01147                 Number current_terms = 
01148                   (this->get_adjoint_rhs(i).dot(this->get_sensitivity_solution(l)) -
01149                    tempvec->dot(this->get_adjoint_solution(i))) / (2.*delta_p);
01150                 sensitivities.second_derivative(i,k,l) += current_terms;
01151 
01152                 // We use the _uk terms twice; symmetry lets us reuse
01153                 // these calculations for the _ul terms.
01154 
01155                 sensitivities.second_derivative(i,l,k) += current_terms;
01156               }
01157         }
01158  
01159       *parameters[k] = old_parameterk - delta_p;
01160       this->assembly(false, true);
01161       this->assemble_qoi_derivative(qoi_indices);
01162 
01163       for (unsigned int l=0; l != Np; ++l)
01164         {
01165           this->matrix->vector_mult(*tempvec, this->get_sensitivity_solution(l));
01166           for (unsigned int i=0; i != Nq; ++i)
01167             if (qoi_indices.has_index(i))
01168               {
01169                 Number current_terms = 
01170                   (-this->get_adjoint_rhs(i).dot(this->get_sensitivity_solution(l)) +
01171                    tempvec->dot(this->get_adjoint_solution(i))) / (2.*delta_p);
01172                 sensitivities.second_derivative(i,k,l) += current_terms;
01173 
01174                 // We use the _uk terms twice; symmetry lets us reuse
01175                 // these calculations for the _ul terms.
01176 
01177                 sensitivities.second_derivative(i,l,k) += current_terms;
01178               }
01179         }
01180 
01181       // Don't leave the parameter perturbed
01182       *parameters[k] = old_parameterk;
01183  
01184       // Our last remaining terms are -R_uu(u,z)*u_k*u_l and
01185       // Q_uu(u)*u_k*u_l
01186       //
01187       // We take directional central finite differences of R_u and Q_u
01188       // to approximate these terms, e.g.:
01189       //
01190       // Q_uu(u)*u_k ~= (Q_u(u+dp*u_k) - Q_u(u-dp*u_k))/(2*dp)
01191 
01192       *this->solution = this->get_sensitivity_solution(k);
01193       *this->solution *= delta_p;
01194       *this->solution += *oldsolution;
01195       this->assembly(false, true);
01196       this->assemble_qoi_derivative(qoi_indices);
01197 
01198       // The Hessian is symmetric, so we just calculate the lower
01199       // triangle and the diagonal, and we get the upper triangle from
01200       // the transpose of the lower
01201       //
01202       // Note that, because we took the directional finite difference
01203       // with respect to k and not l, we've added an O(delta_p^2)
01204       // error to any permutational symmetry in the Hessian...
01205       for (unsigned int l=0; l != k+1; ++l)
01206         {
01207           this->matrix->vector_mult(*tempvec, this->get_sensitivity_solution(l));
01208           for (unsigned int i=0; i != Nq; ++i)
01209             if (qoi_indices.has_index(i))
01210               {
01211                 Number current_terms =
01212                   (this->get_adjoint_rhs(i).dot(this->get_sensitivity_solution(l)) -
01213                    tempvec->dot(this->get_adjoint_solution(i))) / (2.*delta_p);
01214                 sensitivities.second_derivative(i,k,l) += current_terms;
01215                 if (k != l)
01216                   sensitivities.second_derivative(i,l,k) += current_terms;
01217               }
01218         }
01219 
01220       *this->solution = this->get_sensitivity_solution(k);
01221       *this->solution *= -delta_p;
01222       *this->solution += *oldsolution;
01223       this->assembly(false, true);
01224       this->assemble_qoi_derivative(qoi_indices);
01225 
01226       for (unsigned int l=0; l != k+1; ++l)
01227         {
01228           this->matrix->vector_mult(*tempvec, this->get_sensitivity_solution(l));
01229           for (unsigned int i=0; i != Nq; ++i)
01230             if (qoi_indices.has_index(i))
01231               {
01232                 Number current_terms =
01233                   (-this->get_adjoint_rhs(i).dot(this->get_sensitivity_solution(l)) +
01234                    tempvec->dot(this->get_adjoint_solution(i))) / (2.*delta_p);
01235                 sensitivities.second_derivative(i,k,l) += current_terms;
01236                 if (k != l)
01237                   sensitivities.second_derivative(i,l,k) += current_terms;
01238               }
01239         }
01240 
01241       // Don't leave the solution perturbed
01242       *this->solution = *oldsolution;
01243     }
01244 
01245   // All parameters have been reset.
01246   // Don't leave the qoi or system changed - principle of least
01247   // surprise.
01248   this->assembly(true, true);
01249   this->assemble_qoi(qoi_indices);
01250 }

void ImplicitSystem::qoi_parameter_hessian_vector_product	(	const QoISet &	qoi_indices,
		const ParameterVector &	parameters,
		const ParameterVector &	vector,
		SensitivityData &	product
	)			`[virtual, inherited]`

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

Reimplemented from System.

Definition at line 825 of file implicit_system.C.

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

00829 {
00830   // We currently get partial derivatives via finite differencing
00831   const Real delta_p = TOLERANCE;
00832 
00833   // We'll use a single temporary vector for matrix-vector-vector products
00834   AutoPtr<NumericVector<Number> > tempvec = this->solution->zero_clone();
00835 
00836   const unsigned int Np = parameters.size();
00837   const unsigned int Nq = qoi.size();
00838 
00839   // For each quantity of interest q, the parameter sensitivity
00840   // Hessian is defined as q''_{kl} = {d^2 q}/{d p_k d p_l}.
00841   // Given a vector of parameter perturbation weights w_l, this
00842   // function evaluates the hessian-vector product sum_l(q''_{kl}*w_l)
00843   //
00844   // We calculate it from values and partial derivatives of the
00845   // quantity of interest function Q, solution u, adjoint solution z,
00846   // parameter sensitivity adjoint solutions z^l, and residual R, as:
00847   //
00848   // sum_l(q''_{kl}*w_l) =
00849   // sum_l(w_l * Q''_{kl}) + Q''_{uk}(u)*(sum_l(w_l u'_l)) -
00850   // R'_k(u, sum_l(w_l*z^l)) - R'_{uk}(u,z)*(sum_l(w_l u'_l) -
00851   // sum_l(w_l*R''_{kl}(u,z))
00852   //
00853   // See the adjoints model document for more details.
00854 
00855   // We first do an adjoint solve to get z for each quantity of
00856   // interest
00857 
00858   this->adjoint_solve(qoi_indices);
00859 
00860   // Get ready to fill in senstivities:
00861   sensitivities.allocate_data(qoi_indices, *this, parameters);
00862 
00863   // We can't solve for all the solution sensitivities u'_l or for all
00864   // of the parameter sensitivity adjoint solutions z^l without
00865   // requiring O(Nq*Np) linear solves.  So we'll solve directly for their
00866   // weighted sum - this is just O(Nq) solves.
00867 
00868   // First solve for sum_l(w_l u'_l).
00869   this->weighted_sensitivity_solve(parameters, vector);
00870 
00871   // Then solve for sum_l(w_l z^l).
00872   this->weighted_sensitivity_adjoint_solve(parameters, vector, qoi_indices);
00873 
00874   for (unsigned int k=0; k != Np; ++k)
00875     {
00876       // We approximate sum_l(w_l * Q''_{kl}) with a central
00877       // differencing pertubation:
00878       // sum_l(w_l * Q''_{kl}) ~= 
00879       // (Q(p + dp*w_l*e_l + dp*e_k) - Q(p - dp*w_l*e_l + dp*e_k) -
00880       // Q(p + dp*w_l*e_l - dp*e_k) + Q(p - dp*w_l*e_l - dp*e_k))/(4*dp^2)
00881 
00882       // The sum(w_l*R''_kl) term requires the same sort of pertubation,
00883       // and so we subtract it in at the same time:
00884       // sum_l(w_l * R''_{kl}) ~= 
00885       // (R(p + dp*w_l*e_l + dp*e_k) - R(p - dp*w_l*e_l + dp*e_k) -
00886       // R(p + dp*w_l*e_l - dp*e_k) + R(p - dp*w_l*e_l - dp*e_k))/(4*dp^2)
00887 
00888       ParameterVector oldparameters, parameterperturbation;
00889       parameters.deep_copy(oldparameters);
00890       vector.deep_copy(parameterperturbation);
00891       parameterperturbation *= delta_p;
00892       parameters += parameterperturbation;
00893 
00894       Number old_parameter = *parameters[k];
00895 
00896       *parameters[k] = old_parameter + delta_p;
00897       this->assemble_qoi(qoi_indices);
00898       this->assembly(true, false);
00899       std::vector<Number> partial2q_term = this->qoi;
00900       std::vector<Number> partial2R_term(this->qoi.size());
00901       for (unsigned int i=0; i != Nq; ++i)
00902         if (qoi_indices.has_index(i))
00903           partial2R_term[i] = this->rhs->dot(this->get_adjoint_solution(i));
00904 
00905       *parameters[k] = old_parameter - delta_p;
00906       this->assemble_qoi(qoi_indices);
00907       this->assembly(true, false);
00908       for (unsigned int i=0; i != Nq; ++i)
00909         if (qoi_indices.has_index(i))
00910           {
00911             partial2q_term[i] += this->qoi[i];
00912             partial2R_term[i] += this->rhs->dot(this->get_adjoint_solution(i));
00913           }
00914 
00915       oldparameters.value_copy(parameters);
00916       parameterperturbation *= -1.0;
00917       parameters += parameterperturbation;
00918 
00919       // Re-center old_parameter, which may be affected by vector
00920       old_parameter = *parameters[k];
00921 
00922       *parameters[k] = old_parameter + delta_p;
00923       this->assemble_qoi(qoi_indices);
00924       this->assembly(true, false);
00925       for (unsigned int i=0; i != Nq; ++i)
00926         if (qoi_indices.has_index(i))
00927           {
00928             partial2q_term[i] += this->qoi[i];
00929             partial2R_term[i] += this->rhs->dot(this->get_adjoint_solution(i));
00930           }
00931 
00932       *parameters[k] = old_parameter - delta_p;
00933       this->assemble_qoi(qoi_indices);
00934       this->assembly(true, false);
00935       for (unsigned int i=0; i != Nq; ++i)
00936         if (qoi_indices.has_index(i))
00937           {
00938             partial2q_term[i] += this->qoi[i];
00939             partial2R_term[i] += this->rhs->dot(this->get_adjoint_solution(i));
00940           }
00941 
00942       for (unsigned int i=0; i != Nq; ++i)
00943         if (qoi_indices.has_index(i))
00944           {
00945             partial2q_term[i] /= (4. * delta_p * delta_p);
00946             partial2R_term[i] /= (4. * delta_p * delta_p);
00947           }
00948 
00949       for (unsigned int i=0; i != Nq; ++i)
00950         if (qoi_indices.has_index(i))
00951           sensitivities[i][k] = partial2q_term[i] - partial2R_term[i];
00952 
00953       // We get (partial q / partial u), R, and 
00954       // (partial R / partial u) from the user, but centrally
00955       // difference to get q_uk, R_k, and R_uk terms:
00956       // (partial R / partial k)
00957       // 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)
00958       // (partial^2 q / partial u partial k)
00959       // q_uk = (q_u(p+dp*e_k) - q_u(p-dp*e_k))/(2*dp)
00960       // (partial^2 R / partial u partial k)
00961       // 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)
00962 
00963       // To avoid creating Nq temporary vectors for q_uk or R_uk, we add
00964       // subterms to the sensitivities output one by one.
00965       //
00966       // FIXME: this is probably a bad order of operations for
00967       // controlling floating point error.
00968 
00969       *parameters[k] = old_parameter + delta_p;
00970       this->assembly(true, true);
00971       this->assemble_qoi_derivative(qoi_indices);
00972 
00973       this->matrix->vector_mult(*tempvec, this->get_weighted_sensitivity_solution());
00974 
00975       for (unsigned int i=0; i != Nq; ++i)
00976         if (qoi_indices.has_index(i))
00977           sensitivities[i][k] += (this->get_adjoint_rhs(i).dot(this->get_weighted_sensitivity_solution()) -
00978                                   this->rhs->dot(this->get_weighted_sensitivity_adjoint_solution(i)) -
00979                                   this->get_adjoint_solution(i).dot(*tempvec)) / (2.*delta_p);
00980  
00981       *parameters[k] = old_parameter - delta_p;
00982       this->assembly(true, true);
00983       this->assemble_qoi_derivative(qoi_indices);
00984 
00985       this->matrix->vector_mult(*tempvec, this->get_weighted_sensitivity_solution());
00986 
00987       for (unsigned int i=0; i != Nq; ++i)
00988         if (qoi_indices.has_index(i))
00989           sensitivities[i][k] += (-this->get_adjoint_rhs(i).dot(this->get_weighted_sensitivity_solution()) +
00990                                   this->rhs->dot(this->get_weighted_sensitivity_adjoint_solution(i)) +
00991                                   this->get_adjoint_solution(i).dot(*tempvec)) / (2.*delta_p);
00992     }
00993 
00994   // All parameters have been reset.
00995   // Don't leave the qoi or system changed - principle of least
00996   // surprise.
00997   this->assembly(true, true);
00998   this->assemble_qoi(qoi_indices);
00999 }

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

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

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

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

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

Definition at line 383 of file system.C.

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

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

void System::re_update ( ) [virtual, inherited]

Re-update the local values when the mesh has changed. This method takes the data updated by update() and makes it up-to-date on the current mesh.

Definition at line 314 of file system.C.

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

00315 {
00316   //const std::vector<unsigned int>& send_list = _dof_map->get_send_list ();
00317 
00318   // If this system is empty... don't do anything!
00319   if(!this->n_vars())
00320     return;
00321 
00322   // Explicitly build a send_list
00323   std::vector<unsigned int> send_list(solution->size());
00324   Utility::iota (send_list.begin(), send_list.end(), 0);
00325   
00326   // Check sizes
00327   libmesh_assert (current_local_solution->size()       == solution->size());
00328   // Not true with ghosted vectors
00329   // libmesh_assert (current_local_solution->local_size() == solution->size());
00330   libmesh_assert (!send_list.empty());
00331   libmesh_assert (send_list.size() <= solution->size());
00332 
00333   // Create current_local_solution from solution.  This will
00334   // put a local copy of solution into current_local_solution.
00335   solution->localize (*current_local_solution, send_list);
00336 }

void System::read_header	(	Xdr &	io,
		const std::string &	version,
		const bool	read_header = `true`,
		const bool	read_additional_data = `true`,
		const bool	read_legacy_format = `false`
	)			`[inherited]`

Reads the basic data header for this System.

Definition at line 78 of file system_io.C.

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

Referenced by EquationSystems::_read_impl().

00083 {
00084   // This method implements the input of a
00085   // System object, embedded in the output of
00086   // an EquationSystems<T_sys>.  This warrants some 
00087   // documentation.  The output file essentially
00088   // consists of 5 sections:
00089   //
00090   // for this system
00091   //
00092   //   5.) The number of variables in the system (unsigned int)
00093   //
00094   //   for each variable in the system
00095   //     
00096   //     6.) The name of the variable (string)
00097   //     
00098   //     7.) Combined in an FEType:
00099   //         - The approximation order(s) of the variable 
00100   //           (Order Enum, cast to int/s)
00101   //         - The finite element family/ies of the variable 
00102   //           (FEFamily Enum, cast to int/s)
00103   // 
00104   //   end variable loop
00105   //
00106   //   8.) The number of additional vectors (unsigned int),      
00107   //
00108   //     for each additional vector in the system object
00109   // 
00110   //     9.) the name of the additional vector  (string)
00111   //
00112   // end system
00113   libmesh_assert (io.reading());
00114   
00115   // Possibly clear data structures and start from scratch.
00116   if (read_header)
00117     this->clear ();
00118 
00119   // Figure out if we need to read infinite element information.
00120   // This will be true if the version string contains " with infinite elements"
00121   const bool read_ifem_info =
00122     (version.rfind(" with infinite elements") < version.size()) ||
00123     libMesh::on_command_line ("--read_ifem_systems");
00124   
00125   
00126   {
00127     // 5.) 
00128     // Read the number of variables in the system
00129     unsigned int n_vars=0;
00130     if (libMesh::processor_id() == 0) io.data (n_vars);
00131     Parallel::broadcast(n_vars);
00132       
00133     for (unsigned int var=0; var<n_vars; var++)
00134       {             
00135         // 6.)
00136         // Read the name of the var-th variable
00137         std::string var_name;     
00138         if (libMesh::processor_id() == 0) io.data (var_name);
00139         Parallel::broadcast(var_name);
00140                 
00141         // 7.)
00142         // Read the approximation order(s) of the var-th variable 
00143         int order=0;      
00144         if (libMesh::processor_id() == 0) io.data (order);
00145         Parallel::broadcast(order);
00146         
00147         
00148         // do the same for infinite element radial_order 
00149         int rad_order=0;
00150         if (read_ifem_info)
00151           {
00152             if (libMesh::processor_id() == 0) io.data(rad_order);
00153             Parallel::broadcast(rad_order);
00154           }
00155 
00156 
00157         // Read the finite element type of the var-th variable 
00158         int fam=0;
00159         if (libMesh::processor_id() == 0) io.data (fam);
00160         Parallel::broadcast(fam);
00161         FEType type;
00162         type.order  = static_cast<Order>(order);
00163         type.family = static_cast<FEFamily>(fam);
00164 
00165         // Check for incompatibilities.  The shape function indexing was
00166         // changed for the monomial and xyz finite element families to 
00167         // simplify extension to arbitrary p.  The consequence is that 
00168         // old restart files will not be read correctly.  This is expected
00169         // to be an unlikely occurance, but catch it anyway.
00170         if (read_legacy_format)
00171           if ((type.family == MONOMIAL || type.family == XYZ) && 
00172               ((type.order > 2 && this->get_mesh().mesh_dimension() == 2) ||
00173                (type.order > 1 && this->get_mesh().mesh_dimension() == 3)))
00174             {
00175               libmesh_here();
00176               std::cout << "*****************************************************************\n"
00177                         << "* WARNING: reading a potentially incompatible restart file!!!   *\n"
00178                         << "*  contact libmesh-users@lists.sourceforge.net for more details *\n"
00179                         << "*****************************************************************"
00180                         << std::endl;
00181             }
00182                         
00183         // Read additional information for infinite elements    
00184         int radial_fam=0;
00185         int i_map=0;
00186         if (read_ifem_info)
00187           {
00188             if (libMesh::processor_id() == 0) io.data (radial_fam);
00189             Parallel::broadcast(radial_fam);
00190             if (libMesh::processor_id() == 0) io.data (i_map);
00191             Parallel::broadcast(i_map);
00192           }
00193         
00194 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
00195         
00196         type.radial_order  = static_cast<Order>(rad_order);
00197         type.radial_family = static_cast<FEFamily>(radial_fam);
00198         type.inf_map       = static_cast<InfMapType>(i_map);      
00199 
00200 #endif
00201 
00202         if (read_header) 
00203           this->add_variable (var_name, type);
00204       }
00205   }
00206 
00207   // 8.)  
00208   // Read the number of additional vectors.  
00209   unsigned int n_vectors=0;  
00210   if (libMesh::processor_id() == 0) io.data (n_vectors);
00211   Parallel::broadcast(n_vectors);
00212   
00213   // If n_vectors > 0, this means that write_additional_data
00214   // was true when this file was written.  We will need to
00215   // make use of this fact later.
00216   if (n_vectors > 0)
00217     this->_additional_data_written = true;  
00218   
00219   for (unsigned int vec=0; vec<n_vectors; vec++)
00220     {
00221       // 9.)
00222       // Read the name of the vec-th additional vector
00223       std::string vec_name;      
00224       if (libMesh::processor_id() == 0) io.data (vec_name);
00225       Parallel::broadcast(vec_name);
00226       
00227       if (read_additional_data)
00228         {
00229           // sanity checks
00230           libmesh_assert(this->_can_add_vectors);
00231           // Some systems may have added their own vectors already
00232 //        libmesh_assert(this->_vectors.count(vec_name) == 0);
00233 
00234           this->add_vector(vec_name);
00235         }
00236     }
00237 }

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

Reads additional data, namely vectors, for this System.

Definition at line 241 of file system_io.C.

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

00243 {
00244   libmesh_deprecated();
00245   
00246   // This method implements the output of the vectors
00247   // contained in this System object, embedded in the 
00248   // output of an EquationSystems<T_sys>. 
00249   //
00250   //   10.) The global solution vector, re-ordered to be node-major 
00251   //       (More on this later.)                                    
00252   //                                                                
00253   //      for each additional vector in the object          
00254   //                                                                
00255   //      11.) The global additional vector, re-ordered to be       
00256   //           node-major (More on this later.)
00257   libmesh_assert (io.reading());
00258 
00259   // read and reordering buffers
00260   std::vector<Number> global_vector;
00261   std::vector<Number> reordered_vector;
00262 
00263   // 10.)
00264   // Read and set the solution vector
00265   {     
00266     if (libMesh::processor_id() == 0) io.data (global_vector);    
00267     Parallel::broadcast(global_vector);
00268     
00269     // Remember that the stored vector is node-major.
00270     // We need to put it into whatever application-specific
00271     // ordering we may have using the dof_map.
00272     reordered_vector.resize(global_vector.size());
00273 
00274     //std::cout << "global_vector.size()=" << global_vector.size() << std::endl;
00275     //std::cout << "this->n_dofs()=" << this->n_dofs() << std::endl;
00276     
00277     libmesh_assert (global_vector.size() == this->n_dofs());
00278         
00279     unsigned int cnt=0;
00280 
00281     const unsigned int sys     = this->number();
00282     const unsigned int n_vars  = this->n_vars();
00283 
00284     for (unsigned int var=0; var<n_vars; var++)
00285       { 
00286         // First reorder the nodal DOF values
00287         {
00288           MeshBase::node_iterator
00289             it  = this->get_mesh().nodes_begin(),
00290             end = this->get_mesh().nodes_end();
00291         
00292           for (; it != end; ++it)
00293             for (unsigned int index=0; index<(*it)->n_comp(sys,var); index++)
00294               {
00295                 libmesh_assert ((*it)->dof_number(sys, var, index) !=
00296                         DofObject::invalid_id);
00297                 
00298                 libmesh_assert (cnt < global_vector.size());
00299                 
00300                 reordered_vector[(*it)->dof_number(sys, var, index)] =
00301                   global_vector[cnt++]; 
00302               }
00303         }
00304         
00305         // Then reorder the element DOF values
00306         {
00307           MeshBase::element_iterator
00308             it  = this->get_mesh().active_elements_begin(),
00309             end = this->get_mesh().active_elements_end();
00310 
00311           for (; it != end; ++it)
00312             for (unsigned int index=0; index<(*it)->n_comp(sys,var); index++)
00313               {  
00314                 libmesh_assert ((*it)->dof_number(sys, var, index) !=
00315                         DofObject::invalid_id);
00316                 
00317                 libmesh_assert (cnt < global_vector.size());
00318                 
00319                 reordered_vector[(*it)->dof_number(sys, var, index)] =
00320                   global_vector[cnt++]; 
00321               }
00322         }
00323       }
00324             
00325     *(this->solution) = reordered_vector;
00326   }
00327 
00328   // For each additional vector, simply go through the list.
00329   // ONLY attempt to do this IF additional data was actually
00330   // written to the file for this system (controlled by the
00331   // _additional_data_written flag).  
00332   if (this->_additional_data_written)
00333     {
00334       std::map<std::string, NumericVector<Number>* >::iterator
00335         pos = this->_vectors.begin();
00336   
00337       for (; pos != this->_vectors.end(); ++pos)
00338         {
00339           // 11.)          
00340           // Read the values of the vec-th additional vector.
00341           // Prior do _not_ clear, but fill with zero, since the
00342           // additional vectors _have_ to have the same size
00343           // as the solution vector
00344           std::fill (global_vector.begin(), global_vector.end(), libMesh::zero);
00345 
00346           if (libMesh::processor_id() == 0) io.data (global_vector);
00347           Parallel::broadcast(global_vector);
00348 
00349           // If read_additional_data==true, then we will keep this vector, otherwise
00350           // we are going to throw it away.
00351           if (read_additional_data)
00352             {
00353               // Remember that the stored vector is node-major.
00354               // We need to put it into whatever application-specific
00355               // ordering we may have using the dof_map.
00356               std::fill (reordered_vector.begin(),
00357                          reordered_vector.end(),
00358                          libMesh::zero);
00359         
00360               reordered_vector.resize(global_vector.size());    
00361 
00362               libmesh_assert (global_vector.size() == this->n_dofs());
00363         
00364               unsigned int cnt=0;
00365 
00366               const unsigned int sys     = this->number();
00367               const unsigned int n_vars  = this->n_vars();
00368         
00369               for (unsigned int var=0; var<n_vars; var++)
00370                 {
00371                   // First reorder the nodal DOF values
00372                   {
00373                     MeshBase::node_iterator
00374                       it  = this->get_mesh().nodes_begin(),
00375                       end = this->get_mesh().nodes_end();
00376 
00377                     for (; it!=end; ++it)
00378                       for (unsigned int index=0; index<(*it)->n_comp(sys,var); index++)
00379                         {
00380                           libmesh_assert ((*it)->dof_number(sys, var, index) !=
00381                                   DofObject::invalid_id);
00382                           
00383                           libmesh_assert (cnt < global_vector.size());
00384                           
00385                           reordered_vector[(*it)->dof_number(sys, var, index)] =
00386                             global_vector[cnt++]; 
00387                         }
00388                   }
00389 
00390                   // Then reorder the element DOF values
00391                   {
00392                     MeshBase::element_iterator
00393                       it  = this->get_mesh().active_elements_begin(),
00394                       end = this->get_mesh().active_elements_end();
00395 
00396                     for (; it!=end; ++it)
00397                       for (unsigned int index=0; index<(*it)->n_comp(sys,var); index++)
00398                         {  
00399                           libmesh_assert ((*it)->dof_number(sys, var, index) !=
00400                                   DofObject::invalid_id);
00401                           
00402                           libmesh_assert (cnt < global_vector.size());
00403                           
00404                           reordered_vector[(*it)->dof_number(sys, var, index)] =
00405                             global_vector[cnt++]; 
00406                         }
00407                   }
00408                 }
00409             
00410               // use the overloaded operator=(std::vector) to assign the values
00411               *(pos->second) = reordered_vector;
00412             }
00413         }
00414     } // end if (_additional_data_written)    
00415 }

void System::read_parallel_data	(	Xdr &	io,
		const bool	read_additional_data
	)			`[inherited]`

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

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

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

for each additional vector in the object

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

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

Definition at line 419 of file system_io.C.

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

00421 {
00441   libmesh_assert (io.reading());
00442   libmesh_assert (io.is_open());
00443 
00444   // build the ordered nodes and element maps.
00445   // when writing/reading parallel files we need to iterate
00446   // over our nodes/elements in order of increasing global id().
00447   // however, this is not guaranteed to be ordering we obtain
00448   // by using the node_iterators/element_iterators directly.
00449   // so build a set, sorted by id(), that provides the ordering.
00450   // further, for memory economy build the set but then transfer
00451   // its contents to vectors, which will be sorted.
00452   std::vector<const DofObject*> ordered_nodes, ordered_elements;
00453   {    
00454     std::set<const DofObject*, CompareDofObjectsByID>
00455       ordered_nodes_set (this->get_mesh().local_nodes_begin(),
00456                          this->get_mesh().local_nodes_end());
00457       
00458       ordered_nodes.insert(ordered_nodes.end(),
00459                            ordered_nodes_set.begin(),
00460                            ordered_nodes_set.end());
00461   }
00462   {
00463     std::set<const DofObject*, CompareDofObjectsByID>
00464       ordered_elements_set (this->get_mesh().local_elements_begin(),
00465                             this->get_mesh().local_elements_end());
00466       
00467       ordered_elements.insert(ordered_elements.end(),
00468                               ordered_elements_set.begin(),
00469                               ordered_elements_set.end());
00470   }
00471   
00472   std::vector<Number> io_buffer;
00473   
00474   // 9.)
00475   //
00476   // Actually read the solution components
00477   // for the ith system to disk
00478   io.data(io_buffer);
00479   
00480   const unsigned int sys_num = this->number();
00481   const unsigned int n_vars  = this->n_vars();
00482 
00483   unsigned int cnt=0;
00484   
00485   // Loop over each variable and each node, and read out the value.
00486   for (unsigned int var=0; var<n_vars; var++)
00487     {
00488       // First read the node DOF values
00489       for (std::vector<const DofObject*>::const_iterator
00490              it = ordered_nodes.begin(); it != ordered_nodes.end(); ++it)
00491         for (unsigned int comp=0; comp<(*it)->n_comp(sys_num, var); comp++)
00492           {
00493             libmesh_assert ((*it)->dof_number(sys_num, var, comp) !=
00494                     DofObject::invalid_id);
00495             libmesh_assert (cnt < io_buffer.size());
00496             this->solution->set((*it)->dof_number(sys_num, var, comp), io_buffer[cnt++]);
00497           }
00498 
00499       // Then read the element DOF values
00500       for (std::vector<const DofObject*>::const_iterator
00501              it = ordered_elements.begin(); it != ordered_elements.end(); ++it)
00502         for (unsigned int comp=0; comp<(*it)->n_comp(sys_num, var); comp++)
00503           {
00504             libmesh_assert ((*it)->dof_number(sys_num, var, comp) !=
00505                     DofObject::invalid_id);
00506             libmesh_assert (cnt < io_buffer.size());
00507             this->solution->set((*it)->dof_number(sys_num, var, comp), io_buffer[cnt++]);
00508           }      
00509     }
00510   
00511   // Only read additional vectors if wanted  
00512   if (read_additional_data)
00513     {     
00514       std::map<std::string, NumericVector<Number>* >::const_iterator
00515         pos = _vectors.begin();
00516   
00517       for(; pos != this->_vectors.end(); ++pos)
00518         {
00519           cnt=0;
00520           io_buffer.clear();
00521           
00522           // 10.)
00523           //
00524           // Actually read the additional vector components
00525           // for the ith system to disk
00526           io.data(io_buffer);
00527           
00528           // Loop over each variable and each node, and read out the value.
00529           for (unsigned int var=0; var<n_vars; var++)
00530             {
00531               // First read the node DOF values
00532               for (std::vector<const DofObject*>::const_iterator
00533                      it = ordered_nodes.begin(); it != ordered_nodes.end(); ++it)
00534                 for (unsigned int comp=0; comp<(*it)->n_comp(sys_num, var); comp++)
00535                   {
00536                     libmesh_assert ((*it)->dof_number(sys_num, var, comp) !=
00537                             DofObject::invalid_id);
00538                     libmesh_assert (cnt < io_buffer.size());
00539                     this->solution->set((*it)->dof_number(sys_num, var, comp), io_buffer[cnt++]);
00540                   }          
00541               
00542               // Then read the element DOF values
00543               for (std::vector<const DofObject*>::const_iterator
00544                      it = ordered_elements.begin(); it != ordered_elements.end(); ++it)
00545                 for (unsigned int comp=0; comp<(*it)->n_comp(sys_num, var); comp++)
00546                   {
00547                     libmesh_assert ((*it)->dof_number(sys_num, var, comp) !=
00548                             DofObject::invalid_id);
00549                     libmesh_assert (cnt < io_buffer.size());
00550                     this->solution->set((*it)->dof_number(sys_num, var, comp), io_buffer[cnt++]);
00551                   }           
00552             }
00553         }
00554     }
00555 }

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

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

Definition at line 559 of file system_io.C.

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

00561 {
00562   // This method implements the input of the vectors
00563   // contained in this System object, embedded in the 
00564   // output of an EquationSystems<T_sys>. 
00565   //
00566   //   10.) The global solution vector, re-ordered to be node-major 
00567   //       (More on this later.)                                    
00568   //                                                                
00569   //      for each additional vector in the object          
00570   //                                                                
00571   //      11.) The global additional vector, re-ordered to be       
00572   //          node-major (More on this later.)
00573   parallel_only();
00574   std::string comment;
00575 
00576   // 10.)
00577   // Read the global solution vector
00578   {
00579     this->read_serialized_vector(io, *this->solution); 
00580 
00581     // get the comment
00582     if (libMesh::processor_id() == 0)
00583       io.comment (comment);  
00584   }
00585   
00586   // 11.)
00587   // Only read additional vectors if wanted  
00588   if (read_additional_data)
00589     {     
00590       std::map<std::string, NumericVector<Number>* >::const_iterator
00591         pos = _vectors.begin();
00592   
00593       for(; pos != this->_vectors.end(); ++pos)
00594         {
00595           this->read_serialized_vector(io, *pos->second);
00596 
00597           // get the comment
00598           if (libMesh::processor_id() == 0)
00599             io.comment (comment);         
00600             
00601         }
00602     }
00603 }

void DifferentiableSystem::reinit ( ) [virtual, inherited]

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

Reimplemented from ImplicitSystem.

Definition at line 71 of file diff_system.C.

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

00072 {
00073   Parent::reinit();
00074 
00075   time_solver->reinit();
00076 }

void DifferentiableSystem::release_linear_solver ( LinearSolver< Number > * s ) const [virtual, inherited]

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

Reimplemented from ImplicitSystem.

Definition at line 131 of file diff_system.C.

00132 {
00133 }

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

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

Definition at line 213 of file implicit_system.C.

References ImplicitSystem::_matrices.

00214 {
00215   // Make sure the matrix exists
00216   matrices_iterator pos = _matrices.find (mat_name);
00217   
00218   if (pos == _matrices.end())
00219     return NULL;
00220   
00221   return pos->second;
00222 }

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

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

Definition at line 200 of file implicit_system.C.

References ImplicitSystem::_matrices.

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

00201 {
00202   // Make sure the matrix exists
00203   const_matrices_iterator pos = _matrices.find (mat_name);
00204   
00205   if (pos == _matrices.end())
00206     return NULL;
00207   
00208   return pos->second;
00209 }

NumericVector< Number > * System::request_vector ( const unsigned int vec_num ) [inherited]

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

Definition at line 611 of file system.C.

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

00612 {
00613   vectors_iterator v = vectors_begin();
00614   vectors_iterator v_end = vectors_end();
00615   unsigned int num = 0;
00616   while((num<vec_num) && (v!=v_end))
00617     {
00618       num++;
00619       ++v;
00620     }
00621   if (v==v_end)
00622     return NULL;
00623   return v->second;
00624 }

const NumericVector< Number > * System::request_vector ( const unsigned int vec_num ) const [inherited]

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

Definition at line 594 of file system.C.

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

00595 {
00596   const_vectors_iterator v = vectors_begin();
00597   const_vectors_iterator v_end = vectors_end();
00598   unsigned int num = 0;
00599   while((num<vec_num) && (v!=v_end))
00600     {
00601       num++;
00602       ++v;
00603     }
00604   if (v==v_end)
00605     return NULL;
00606   return v->second;
00607 }

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

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

Definition at line 582 of file system.C.

References System::_vectors.

00583 {
00584   vectors_iterator pos = _vectors.find(vec_name);
00585   
00586   if (pos == _vectors.end())
00587     return NULL;
00588   
00589   return pos->second;
00590 }

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

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

Definition at line 570 of file system.C.

References System::_vectors.

00571 {
00572   const_vectors_iterator pos = _vectors.find(vec_name);
00573   
00574   if (pos == _vectors.end())
00575     return NULL;
00576   
00577   return pos->second;
00578 }

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

Restrict vectors after the mesh has coarsened

Definition at line 219 of file system.C.

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

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

00220 {
00221 #ifdef LIBMESH_ENABLE_AMR
00222   // Restrict the _vectors on the coarsened cells
00223   for (vectors_iterator pos = _vectors.begin(); pos != _vectors.end(); ++pos)
00224     {
00225       NumericVector<Number>* v = pos->second;
00226       
00227       if (_vector_projections[pos->first])
00228         this->project_vector (*v);
00229       else
00230         v->init (this->n_dofs(), this->n_local_dofs(), false, PARALLEL);
00231     }
00232 
00233   const std::vector<unsigned int>& send_list = _dof_map->get_send_list ();
00234 
00235   // Restrict the solution on the coarsened cells
00236   if (_solution_projection)
00237     this->project_vector (*solution);
00238 
00239 #ifdef LIBMESH_ENABLE_GHOSTED
00240   current_local_solution->init(this->n_dofs(),
00241                                this->n_local_dofs(), send_list, 
00242                                false, GHOSTED);
00243 #else
00244   current_local_solution->init(this->n_dofs());
00245 #endif
00246 
00247   if (_solution_projection)
00248     solution->localize (*current_local_solution, send_list); 
00249 
00250 #endif // LIBMESH_ENABLE_AMR
00251 }

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

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

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

Reimplemented from System.

Definition at line 282 of file implicit_system.C.

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

00283 {
00284   // Log how long the linear solve takes.
00285   START_LOG("sensitivity_solve()", "ImplicitSystem");
00286 
00287   // The forward system should now already be solved.
00288   // Now assemble the corresponding sensitivity system.
00289 
00290   if (this->assemble_before_solve)
00291     {
00292       // Build the Jacobian
00293       this->assembly(false, true);
00294       this->matrix->close();
00295 
00296       // Reset and build the RHS from the residual derivatives
00297       this->assemble_residual_derivatives(parameters);
00298     }
00299 
00300   // The sensitivity problem is linear
00301   LinearSolver<Number> *linear_solver = this->get_linear_solver();
00302 
00303   // Our iteration counts and residuals will be sums of the individual
00304   // results
00305   std::pair<unsigned int, Real> solver_params =
00306     this->get_linear_solve_parameters();
00307   std::pair<unsigned int, Real> totalrval = std::make_pair(0,0.0);
00308 
00309   // Solve the linear system.
00310   SparseMatrix<Number> *pc = this->request_matrix("Preconditioner");
00311   for (unsigned int p=0; p != parameters.size(); ++p)
00312     {
00313       std::pair<unsigned int, Real> rval =
00314         linear_solver->solve (*matrix, pc,
00315                               this->get_sensitivity_solution(p),
00316                               this->get_sensitivity_rhs(p),
00317                               solver_params.second,
00318                               solver_params.first);
00319 
00320       totalrval.first  += rval.first;
00321       totalrval.second += rval.second;
00322     }
00323 
00324   // The linear solver may not have fit our constraints exactly
00325 #ifdef LIBMESH_ENABLE_AMR
00326   for (unsigned int p=0; p != parameters.size(); ++p)
00327     this->get_dof_map().enforce_constraints_exactly
00328       (*this, &this->get_sensitivity_solution(p));
00329 #endif
00330 
00331   this->release_linear_solver(linear_solver);
00332 
00333   // Stop logging the nonlinear solve
00334   STOP_LOG("sensitivity_solve()", "ImplicitSystem");
00335 
00336   return totalrval;
00337 }

virtual bool DifferentiableSystem::side_constraint	(	bool	request_jacobian,
		DiffContext &
	)			`[inline, virtual, inherited]`

Adds the time derivative contribution on side of elem to elem_residual. If this method receives request_jacobian = true, then it should compute elem_jacobian and return true if possible. If elem_jacobian has not been computed then the method should return false.

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

Definition at line 197 of file diff_system.h.

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

00198                                                {
00199     return request_jacobian;
00200   }

virtual bool DifferentiableSystem::side_mass_residual	(	bool	request_jacobian,
		DiffContext &
	)			`[inline, virtual, inherited]`

Adds a mass vector contribution on side of elem to elem_residual. If this method receives request_jacobian = true, then it should compute elem_jacobian and return true if possible. If elem_jacobian has not been computed then the method should return false.

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

Definition at line 336 of file diff_system.h.

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

00337                                                   {
00338     return request_jacobian;
00339   }

virtual void DifferentiableSystem::side_postprocess ( DiffContext & ) [inline, virtual, inherited]

Does any work that needs to be done on side of elem in a postprocessing loop.

Definition at line 231 of file diff_system.h.

Referenced by postprocess().

00231 {}

virtual void DifferentiableSystem::side_qoi	(	DiffContext &	,
		const QoISet &
	)			`[inline, virtual, inherited]`

Does any work that needs to be done on side of elem in a quantity of interest assembly loop, outputting to elem_qoi.

Only qois included in the supplied QoISet need to be assembled.

Definition at line 263 of file diff_system.h.

Referenced by assemble_qoi().

00265     {}

virtual void DifferentiableSystem::side_qoi_derivative	(	DiffContext &	,
		const QoISet &
	)			`[inline, virtual, inherited]`

Does any work that needs to be done on side of elem in a quantity of interest derivative assembly loop, outputting to elem_qoi_derivative.

Only qois included in the supplied QoISet need their derivatives assembled.

Definition at line 275 of file diff_system.h.

Referenced by assemble_qoi_derivative().

00277     {}

virtual bool DifferentiableSystem::side_time_derivative	(	bool	request_jacobian,
		DiffContext &
	)			`[inline, virtual, inherited]`

Adds the time derivative contribution on side of elem to elem_residual. If this method receives request_jacobian = true, then it should compute elem_jacobian and return true if possible. If elem_jacobian has not been computed then the method should return false.

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

Definition at line 181 of file diff_system.h.

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

00182                                                     {
00183     return request_jacobian;
00184   }

void FEMSystem::solve ( ) [virtual]

Invokes the solver associated with the system. For steady state solvers, this will find a root x where F(x) = 0. For transient solvers, this will integrate dx/dt = F(x).

For moving mesh systems, this also translates the mesh to the solution position.

Reimplemented from DifferentiableSystem.

Reimplemented in ContinuationSystem.

Definition at line 383 of file fem_system.C.

References mesh_position_set(), and DifferentiableSystem::solve().

00384 {
00385   Parent::solve();
00386 
00387   this->mesh_position_set();
00388 }

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

Returns:: a clever pointer to the system.

Reimplemented from ExplicitSystem.

Reimplemented in LinearImplicitSystem, and NonlinearImplicitSystem.

Definition at line 72 of file implicit_system.h.

00072 { return *this; }

virtual std::string ImplicitSystem::system_type ( ) const [inline, virtual, inherited]

Assembles & solves the linear system Ax=b.

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

Reimplemented from ExplicitSystem.

Reimplemented in FrequencySystem, LinearImplicitSystem, NewmarkSystem, and NonlinearImplicitSystem.

Definition at line 107 of file implicit_system.h.

00107 { return "Implicit"; }

void FEMSystem::time_evolving ( unsigned int var ) [virtual]

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

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

Reimplemented from DifferentiableSystem.

Definition at line 766 of file fem_system.C.

References DifferentiableSystem::time_evolving().

00767 {
00768   // Call the parent function
00769   Parent::time_evolving(var);
00770 }

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

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

Definition at line 295 of file system.C.

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

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

00296 {
00297   const std::vector<unsigned int>& send_list = _dof_map->get_send_list ();
00298 
00299   // Check sizes
00300   libmesh_assert (current_local_solution->size() == solution->size());
00301 // More processors than elements => empty send_list
00302 //  libmesh_assert (!send_list.empty());
00303   libmesh_assert (send_list.size() <= solution->size());
00304 
00305   // Create current_local_solution from solution.  This will
00306   // put a local copy of solution into current_local_solution.
00307   // Only the necessary values (specified by the send_list)
00308   // are copied to minimize communication
00309   solution->localize (*current_local_solution, send_list); 
00310 }

void 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 539 of file system.C.

References System::solution.

00541 {
00542   global_soln.resize        (solution->size());
00543 
00544   solution->localize_to_one (global_soln, dest_proc);
00545 }

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

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

Definition at line 530 of file system.C.

References System::solution.

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

00531 {
00532   global_soln.resize (solution->size());
00533 
00534   solution->localize (global_soln);
00535 }

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

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

Definition at line 1367 of file system.C.

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

Referenced by System::assemble().

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

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

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

Definition at line 1377 of file system.C.

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

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

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

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

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

Definition at line 1358 of file system.C.

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

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

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

void 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 1387 of file system.C.

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

Referenced by System::assemble_qoi().

01388 {
01389   // Call the user-provided quantity of interest function, 
01390   // if it was provided
01391   if(_qoi_evaluate != NULL)
01392     this->_qoi_evaluate(_equation_systems, this->name(), qoi_indices);
01393 }

void 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 1397 of file system.C.

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

Referenced by System::assemble_qoi_derivative().

01398 {
01399   // Call the user-provided quantity of interest derivative, 
01400   // if it was provided
01401   if(_qoi_evaluate_derivative != NULL)
01402     this->_qoi_evaluate_derivative(_equation_systems, this->name(), qoi_indices);
01403 }

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

Return a constant reference to Variable var.

Definition at line 1404 of file system.h.

References System::_variables.

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

01405 {
01406   libmesh_assert (i < _variables.size());
01407 
01408   return _variables[i];
01409 }

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

Returns:: the name of variable i.

Definition at line 1414 of file system.h.

References System::_variables.

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

01415 {
01416   libmesh_assert (i < _variables.size());
01417 
01418   return _variables[i].name();
01419 }

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

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

Definition at line 942 of file system.C.

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

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

00943 {
00944   // Make sure the variable exists
00945   std::map<std::string, unsigned short int>::const_iterator
00946     pos = _variable_numbers.find(var);
00947   
00948   if (pos == _variable_numbers.end())
00949     {
00950       std::cerr << "ERROR: variable "
00951                 << var
00952                 << " does not exist in this system!"
00953                 << std::endl;      
00954       libmesh_error();
00955     }
00956   libmesh_assert (_variables[pos->second].name() == var);
00957 
00958   return pos->second;
00959 }

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

Returns:: the finite element type for variable var.

Definition at line 1434 of file system.h.

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

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

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

Returns:: the finite element type variable number i.

Definition at line 1424 of file system.h.

References System::_variables.

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

01425 {
01426   libmesh_assert (i < _variables.size());
01427   
01428   return _variables[i].type();
01429 }

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

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

Definition at line 698 of file system.C.

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

00699 {
00700   const_vectors_iterator v = vectors_begin();
00701   const_vectors_iterator v_end = vectors_end();
00702   unsigned int num = 0;
00703   while((num<vec_num) && (v!=v_end))
00704     {
00705       num++;
00706       ++v;
00707     }
00708   libmesh_assert(v!=v_end);
00709   return v->first;
00710 }

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

Beginning of vectors container

Definition at line 1470 of file system.h.

References System::_vectors.

01471 {
01472   return _vectors.begin();
01473 }

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

Beginning of vectors container

Definition at line 1464 of file system.h.

References System::_vectors.

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

01465 {
01466   return _vectors.begin();
01467 }

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

End of vectors container

Definition at line 1482 of file system.h.

References System::_vectors.

01483 {
01484   return _vectors.end();
01485 }

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

End of vectors container

Definition at line 1476 of file system.h.

References System::_vectors.

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

01477 {
01478   return _vectors.end();
01479 }

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

Definition at line 404 of file implicit_system.C.

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

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

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

Definition at line 540 of file implicit_system.C.

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

00542 {
00543   // Log how long the linear solve takes.
00544   START_LOG("weighted_sensitivity_solve()", "ImplicitSystem");
00545 
00546   // We currently get partial derivatives via central differencing
00547   const Real delta_p = TOLERANCE;
00548 
00549   // The forward system should now already be solved.
00550 
00551   // Now we're assembling a weighted sum of sensitivity systems:
00552   //
00553   // dR/du (u, v)(sum(w_l*u'_l)) = -sum_l(w_l*R'_l (u, v)) forall v
00554 
00555   // We'll assemble the rhs first, because the R' term will require
00556   // perturbing the system, and some applications may not be able to
00557   // assemble a perturbed residual without simultaneously constructing
00558   // a perturbed jacobian.
00559 
00560   // We approximate the _l partial derivatives via a central
00561   // differencing perturbation in the w_l direction:
00562   //
00563   // sum_l(w_l*v_l) ~= (v(p + dp*w_l*e_l) - v(p - dp*w_l*e_l))/(2*dp)
00564 
00565   ParameterVector oldparameters, parameterperturbation;
00566   parameters.deep_copy(oldparameters);
00567   weights.deep_copy(parameterperturbation);
00568   parameterperturbation *= delta_p;
00569   parameters += parameterperturbation;
00570 
00571   this->assembly(true, false);
00572   this->rhs->close();
00573 
00574   AutoPtr<NumericVector<Number> > temprhs = this->rhs->clone();
00575 
00576   oldparameters.value_copy(parameters);
00577   parameterperturbation *= -1.0;
00578   parameters += parameterperturbation;
00579 
00580   this->assembly(true, false);
00581   this->rhs->close();
00582 
00583   *temprhs -= *(this->rhs);
00584   *temprhs /= (2.0*delta_p);
00585 
00586   // Finally, assemble the jacobian at the non-perturbed parameter
00587   // values
00588   oldparameters.value_copy(parameters);
00589 
00590   // Build the Jacobian
00591   this->assembly(false, true);
00592   this->matrix->close();
00593 
00594   // The weighted sensitivity problem is linear
00595   LinearSolver<Number> *linear_solver = this->get_linear_solver();
00596 
00597   std::pair<unsigned int, Real> solver_params =
00598     this->get_linear_solve_parameters();
00599 
00600   const std::pair<unsigned int, Real> rval =
00601     linear_solver->solve (*matrix, this->add_weighted_sensitivity_solution(),
00602                           *temprhs,
00603                           solver_params.second,
00604                           solver_params.first);
00605 
00606   this->release_linear_solver(linear_solver);
00607 
00608   // The linear solver may not have fit our constraints exactly
00609 #ifdef LIBMESH_ENABLE_AMR
00610   this->get_dof_map().enforce_constraints_exactly
00611     (*this, &this->get_weighted_sensitivity_solution());
00612 #endif
00613 
00614   // Stop logging the nonlinear solve
00615   STOP_LOG("weighted_sensitivity_solve()", "ImplicitSystem");
00616 
00617   return rval;
00618 }

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

7.) Combined in an FEType:

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

end variable loop

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

for each additional vector in the system object

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

end system

Definition at line 815 of file system_io.C.

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

00818 {
00850   libmesh_assert (io.writing());
00851 
00852 
00853   // Only write the header information
00854   // if we are processor 0.
00855   if (this->get_mesh().processor_id() != 0)
00856     return;
00857   
00858   std::string comment;
00859   char buf[80];
00860 
00861   // 5.) 
00862   // Write the number of variables in the system
00863 
00864   {
00865     // set up the comment
00866     comment = "# No. of Variables in System \"";
00867     comment += this->name();
00868     comment += "\"";    
00869           
00870     unsigned int n_vars = this->n_vars();   
00871     io.data (n_vars, comment.c_str());
00872   }
00873 
00874 
00875   for (unsigned int var=0; var<this->n_vars(); var++)
00876     {
00877       // 6.)
00878       // Write the name of the var-th variable
00879       {
00880         // set up the comment
00881         comment  = "#   Name, Variable No. ";
00882         std::sprintf(buf, "%d", var);
00883         comment += buf;
00884         comment += ", System \"";
00885         comment += this->name();
00886         comment += "\"";
00887               
00888         std::string var_name = this->variable_name(var);             
00889         io.data (var_name, comment.c_str());
00890       }
00891       
00892       // 7.)
00893       // Write the approximation order of the var-th variable 
00894       // in this system
00895       {
00896         // set up the comment
00897         comment = "#     Approximation Order, Variable \"";
00898         std::sprintf(buf, "%s", this->variable_name(var).c_str());
00899         comment += buf;
00900         comment += "\", System \"";
00901         comment += this->name();
00902         comment += "\"";
00903               
00904         int order = static_cast<int>(this->variable_type(var).order);         
00905         io.data (order, comment.c_str());
00906       }
00907 
00908 
00909 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
00910       
00911       // do the same for radial_order
00912       {
00913         comment = "#     Radial Approximation Order, Variable \"";
00914         std::sprintf(buf, "%s", this->variable_name(var).c_str());
00915         comment += buf;
00916         comment += "\", System \"";
00917         comment += this->name();
00918         comment += "\"";
00919       
00920         int rad_order = static_cast<int>(this->variable_type(var).radial_order);              
00921         io.data (rad_order, comment.c_str());
00922       }
00923 
00924 #endif
00925       
00926       // Write the Finite Element type of the var-th variable 
00927       // in this System
00928       {
00929         // set up the comment
00930         comment = "#     FE Family, Variable \"";
00931         std::sprintf(buf, "%s", this->variable_name(var).c_str());
00932         comment += buf;
00933         comment += "\", System \"";
00934         comment += this->name();
00935         comment += "\"";
00936       
00937         const FEType& type = this->variable_type(var);        
00938         int fam = static_cast<int>(type.family);              
00939         io.data (fam, comment.c_str());
00940 
00941 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
00942         
00943         comment = "#     Radial FE Family, Variable \"";
00944         std::sprintf(buf, "%s", this->variable_name(var).c_str());
00945         comment += buf;
00946         comment += "\", System \"";
00947         comment += this->name();
00948         comment += "\"";
00949       
00950         int radial_fam = static_cast<int>(type.radial_family);
00951         io.data (radial_fam, comment.c_str());  
00952         
00953         comment = "#     Infinite Mapping Type, Variable \"";
00954         std::sprintf(buf, "%s", this->variable_name(var).c_str());
00955         comment += buf;
00956         comment += "\", System \"";
00957         comment += this->name();
00958         comment += "\"";
00959 
00960         int i_map = static_cast<int>(type.inf_map);           
00961         io.data (i_map, comment.c_str());
00962 #endif
00963       }
00964     } // end of the variable loop
00965  
00966   // 8.) 
00967   // Write the number of additional vectors in the System.
00968   // If write_additional_data==false, then write zero for
00969   // the number of additional vectors.
00970   {       
00971     {
00972       // set up the comment
00973       comment = "# No. of Additional Vectors, System \"";
00974       comment += this->name();
00975       comment += "\"";
00976     
00977       unsigned int n_vectors = write_additional_data ? this->n_vectors () : 0;
00978       io.data (n_vectors, comment.c_str());
00979     }  
00980       
00981     if (write_additional_data)
00982       {
00983         std::map<std::string, NumericVector<Number>* >::const_iterator
00984           vec_pos = this->_vectors.begin();
00985         unsigned int cnt=0;
00986         
00987         for (; vec_pos != this->_vectors.end(); ++vec_pos)
00988           {
00989             // 9.)
00990             // write the name of the cnt-th additional vector
00991             comment =  "# Name of ";
00992             std::sprintf(buf, "%d", cnt++);
00993             comment += buf;
00994             comment += "th vector";
00995             std::string vec_name = vec_pos->first;
00996             
00997             io.data (vec_name, comment.c_str());
00998           }
00999       }
01000   }
01001 }

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

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

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

for each additional vector in the object

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

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

Definition at line 1213 of file system_io.C.

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

01215 {
01235   std::string comment;
01236   
01237   libmesh_assert (io.writing());
01238 
01239   std::vector<Number> io_buffer; io_buffer.reserve(this->solution->local_size());
01240 
01241   // build the ordered nodes and element maps.
01242   // when writing/reading parallel files we need to iterate
01243   // over our nodes/elements in order of increasing global id().
01244   // however, this is not guaranteed to be ordering we obtain
01245   // by using the node_iterators/element_iterators directly.
01246   // so build a set, sorted by id(), that provides the ordering.
01247   // further, for memory economy build the set but then transfer
01248   // its contents to vectors, which will be sorted.
01249   std::vector<const DofObject*> ordered_nodes, ordered_elements;
01250   {    
01251     std::set<const DofObject*, CompareDofObjectsByID>
01252       ordered_nodes_set (this->get_mesh().local_nodes_begin(),
01253                          this->get_mesh().local_nodes_end());
01254       
01255       ordered_nodes.insert(ordered_nodes.end(),
01256                            ordered_nodes_set.begin(),
01257                            ordered_nodes_set.end());
01258   }
01259   {
01260     std::set<const DofObject*, CompareDofObjectsByID>
01261       ordered_elements_set (this->get_mesh().local_elements_begin(),
01262                             this->get_mesh().local_elements_end());
01263       
01264       ordered_elements.insert(ordered_elements.end(),
01265                               ordered_elements_set.begin(),
01266                               ordered_elements_set.end());
01267   }
01268   
01269   const unsigned int sys_num = this->number();
01270   const unsigned int n_vars  = this->n_vars();
01271   
01272   // Loop over each variable and each node, and write out the value.
01273   for (unsigned int var=0; var<n_vars; var++)
01274     {
01275       // First write the node DOF values
01276       for (std::vector<const DofObject*>::const_iterator
01277              it = ordered_nodes.begin(); it != ordered_nodes.end(); ++it)      
01278         for (unsigned int comp=0; comp<(*it)->n_comp(sys_num, var); comp++)
01279           {
01280             //std::cout << "(*it)->id()=" << (*it)->id() << std::endl;
01281             libmesh_assert ((*it)->dof_number(sys_num, var, comp) !=
01282                     DofObject::invalid_id);
01283             
01284             io_buffer.push_back((*this->solution)((*it)->dof_number(sys_num, var, comp)));
01285           }
01286 
01287       // Then write the element DOF values
01288       for (std::vector<const DofObject*>::const_iterator
01289              it = ordered_elements.begin(); it != ordered_elements.end(); ++it)
01290         for (unsigned int comp=0; comp<(*it)->n_comp(sys_num, var); comp++)
01291           {
01292             libmesh_assert ((*it)->dof_number(sys_num, var, comp) !=
01293                     DofObject::invalid_id);
01294               
01295             io_buffer.push_back((*this->solution)((*it)->dof_number(sys_num, var, comp)));    
01296           }
01297     }
01298   
01299   // 9.)
01300   //
01301   // Actually write the reordered solution vector 
01302   // for the ith system to disk
01303   
01304   // set up the comment
01305   {
01306     comment = "# System \"";
01307     comment += this->name();
01308     comment += "\" Solution Vector";
01309   }
01310   
01311   io.data (io_buffer, comment.c_str());   
01312   
01313   // Only write additional vectors if wanted  
01314   if (write_additional_data)
01315     {     
01316       std::map<std::string, NumericVector<Number>* >::const_iterator
01317         pos = _vectors.begin();
01318   
01319       for(; pos != this->_vectors.end(); ++pos)
01320         {
01321           io_buffer.clear(); io_buffer.reserve( pos->second->local_size());
01322           
01323           // Loop over each variable and each node, and write out the value.
01324           for (unsigned int var=0; var<n_vars; var++)
01325             {
01326               // First write the node DOF values
01327               for (std::vector<const DofObject*>::const_iterator
01328                      it = ordered_nodes.begin(); it != ordered_nodes.end(); ++it)      
01329                 for (unsigned int comp=0; comp<(*it)->n_comp(sys_num, var); comp++)
01330                   {
01331                     libmesh_assert ((*it)->dof_number(sys_num, var, comp) !=
01332                             DofObject::invalid_id);
01333                       
01334                     io_buffer.push_back((*pos->second)((*it)->dof_number(sys_num, var, comp)));   
01335                   }
01336 
01337               // Then write the element DOF values
01338               for (std::vector<const DofObject*>::const_iterator
01339                      it = ordered_elements.begin(); it != ordered_elements.end(); ++it)
01340                 for (unsigned int comp=0; comp<(*it)->n_comp(sys_num, var); comp++)
01341                   {
01342                     libmesh_assert ((*it)->dof_number(sys_num, var, comp) !=
01343                             DofObject::invalid_id);
01344               
01345                     io_buffer.push_back((*pos->second)((*it)->dof_number(sys_num, var, comp)));
01346                   }
01347             }
01348           
01349           // 10.)
01350           //
01351           // Actually write the reordered additional vector 
01352           // for this system to disk
01353           
01354           // set up the comment
01355           {
01356             comment = "# System \"";
01357             comment += this->name(); 
01358             comment += "\" Additional Vector \"";
01359             comment += pos->first;
01360             comment += "\"";
01361           }
01362               
01363           io.data (io_buffer, comment.c_str());         
01364         }
01365     }
01366 }

void 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 1370 of file system_io.C.

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

01372 {
01386   parallel_only();
01387   std::string comment;
01388 
01389   this->write_serialized_vector(io, *this->solution); 
01390 
01391   // set up the comment
01392   if (libMesh::processor_id() == 0)
01393     {
01394       comment = "# System \"";
01395       comment += this->name();
01396       comment += "\" Solution Vector";
01397 
01398       io.comment (comment);
01399     }
01400 
01401   // Only write additional vectors if wanted  
01402   if (write_additional_data)
01403     {     
01404       std::map<std::string, NumericVector<Number>* >::const_iterator
01405         pos = _vectors.begin();
01406   
01407       for(; pos != this->_vectors.end(); ++pos)
01408         {
01409           this->write_serialized_vector(io, *pos->second);
01410 
01411           // set up the comment
01412           if (libMesh::processor_id() == 0)
01413             {
01414               comment = "# System \"";
01415               comment += this->name(); 
01416               comment += "\" Additional Vector \"";
01417               comment += pos->first;
01418               comment += "\"";
01419               io.comment (comment);       
01420             }
01421         }
01422     }
01423 }

void 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 995 of file system.C.

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

00996 {
00997   /* Make sure the call makes sense.  */
00998   libmesh_assert(var_num<this->n_vars());
00999 
01000   /* Get a reference to the mesh.  */
01001   const MeshBase& mesh = this->get_mesh();
01002 
01003   /* Check which system we are.  */
01004   const unsigned int sys_num = this->number();
01005 
01006   /* Loop over nodes.  */
01007   {
01008     MeshBase::const_node_iterator it = mesh.local_nodes_begin();
01009     const MeshBase::const_node_iterator end_it = mesh.local_nodes_end(); 
01010     for ( ; it != end_it; ++it)
01011       {    
01012         const Node* node = *it;
01013         unsigned int n_comp = node->n_comp(sys_num,var_num);
01014         for(unsigned int i=0; i<n_comp; i++)
01015           {
01016             const unsigned int index = node->dof_number(sys_num,var_num,i);
01017             v.set(index,0.0);
01018           }
01019       }
01020   }
01021 
01022   /* Loop over elements.  */
01023   {
01024     MeshBase::const_element_iterator it = mesh.active_local_elements_begin();
01025     const MeshBase::const_element_iterator end_it = mesh.active_local_elements_end(); 
01026     for ( ; it != end_it; ++it)
01027       {    
01028         const Elem* elem = *it;
01029         unsigned int n_comp = elem->n_comp(sys_num,var_num);
01030         for(unsigned int i=0; i<n_comp; i++)
01031           {
01032             const unsigned int index = elem->dof_number(sys_num,var_num,i);
01033             v.set(index,0.0);
01034           }
01035       }
01036   }
01037 }

Member Data Documentation

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

Actually holds the data.

Definition at line 110 of file reference_counter.h.

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

unsigned int FEMSystem::_mesh_sys [protected]

System and variables from which to acquire moving mesh information

Definition at line 309 of file fem_system.h.

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

unsigned int FEMSystem::_mesh_x_var [protected]

Definition at line 309 of file fem_system.h.

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

unsigned int FEMSystem::_mesh_y_var [protected]

Definition at line 309 of file fem_system.h.

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

unsigned int FEMSystem::_mesh_z_var [protected]

Definition at line 309 of file fem_system.h.

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

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

Mutual exclusion object to enable thread-safe reference counting.

Definition at line 123 of file reference_counter.h.

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

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

Definition at line 118 of file reference_counter.h.

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

std::vector<bool> DifferentiableSystem::_time_evolving [protected, inherited]

Stores bools to tell us which variables are evolving in time and which are just constraints

Definition at line 443 of file diff_system.h.

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

bool System::assemble_before_solve [inherited]

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

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

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

Definition at line 985 of file system.h.

Referenced by ImplicitSystem::adjoint_solve(), ImplicitSystem::sensitivity_solve(), LinearImplicitSystem::solve(), EigenSystem::solve(), and ImplicitSystem::weighted_sensitivity_adjoint_solve().

bool DifferentiableSystem::assemble_qoi_sides [inherited]

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

Definition at line 220 of file diff_system.h.

Referenced by assemble_qoi(), and assemble_qoi_derivative().

bool DifferentiableSystem::compute_internal_sides [inherited]

compute_internal_sides is false by default, indicating that side_* computations will only be done on boundary sides. If compute_internal_sides is true, computations will be done on sides between elements as well.

Definition at line 168 of file diff_system.h.

Referenced by assemble_qoi(), assemble_qoi_derivative(), assembly(), and postprocess().

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

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

Definition at line 1012 of file system.h.

Referenced by NonlinearImplicitSystem::assembly(), System::clear(), Problem_Interface::computeF(), System::current_solution(), ExactErrorEstimator::estimate_error(), System::init_data(), System::project_solution(), System::re_update(), System::reinit(), System::restrict_vectors(), and System::update().

Real DifferentiableSystem::deltat [inherited]

For time-dependent problems, this is the amount delta t to advance the solution in time.

Definition at line 384 of file diff_system.h.

Referenced by UnsteadySolver::advance_timestep(), EulerSolver::element_residual(), Euler2Solver::element_residual(), eulerian_residual(), EulerSolver::side_residual(), Euler2Solver::side_residual(), UnsteadySolver::solve(), and TwostepTimeSolver::solve().

int FEMSystem::extra_quadrature_order

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 249 of file fem_system.h.

Referenced by FEMContext::FEMContext().

bool FEMSystem::fe_reinit_during_postprocess

If fe_reinit_during_postprocess is true (it is true by default), FE objects will be reinit()ed with their default quadrature rules. If false, FE objects will need to be reinit()ed by the user or will be in an undefined state.

Definition at line 238 of file fem_system.h.

Referenced by postprocess().

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

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

Definition at line 326 of file implicit_system.h.

Referenced by __libmesh_petsc_diff_solver_jacobian(), ImplicitSystem::add_system_matrix(), ImplicitSystem::adjoint_solve(), ImplicitSystem::assemble(), NonlinearImplicitSystem::assembly(), LinearImplicitSystem::assembly(), assembly(), ImplicitSystem::clear(), NewmarkSystem::compute_matrix(), ContinuationSystem::continuation_solve(), ImplicitSystem::init_matrices(), ImplicitSystem::sensitivity_solve(), PetscDiffSolver::solve(), NonlinearImplicitSystem::solve(), NewtonSolver::solve(), LinearImplicitSystem::solve(), FrequencySystem::solve(), EigenTimeSolver::solve(), ContinuationSystem::solve_tangent(), ImplicitSystem::weighted_sensitivity_adjoint_solve(), and ImplicitSystem::weighted_sensitivity_solve().

Real FEMSystem::numerical_jacobian_h

If calculating numeric jacobians is required, the FEMSystem will perturb each solution vector entry by numerical_jacobian_h when calculating finite differences.

Definition at line 256 of file fem_system.h.

Referenced by numerical_jacobian().

bool DifferentiableSystem::postprocess_sides [inherited]

If postprocess_sides is true (it is false by default), the postprocessing loop will loop over all sides as well as all elements.

Definition at line 213 of file diff_system.h.

Referenced by postprocess().

bool DifferentiableSystem::print_element_jacobians [inherited]

Set print_element_jacobians to true to print each J_elem contribution.

Definition at line 421 of file diff_system.h.

Referenced by assembly().

bool DifferentiableSystem::print_jacobian_norms [inherited]

Set print_jacobian_norms to true to print |J| whenever it is assembled.

Definition at line 411 of file diff_system.h.

Referenced by assembly().

bool DifferentiableSystem::print_jacobians [inherited]

Set print_jacobians to true to print J whenever it is assembled.

Definition at line 416 of file diff_system.h.

Referenced by assembly().

bool DifferentiableSystem::print_residual_norms [inherited]

Set print_residual_norms to true to print |F| whenever it is assembled.

Definition at line 401 of file diff_system.h.

Referenced by assembly().

bool DifferentiableSystem::print_residuals [inherited]

Set print_residuals to true to print F whenever it is assembled.

Definition at line 406 of file diff_system.h.

Referenced by assembly().

bool DifferentiableSystem::print_solution_norms [inherited]

Set print_residual_norms to true to print |U| whenever it is used in an assembly() call

Definition at line 390 of file diff_system.h.

Referenced by assembly().

bool DifferentiableSystem::print_solutions [inherited]

Set print_solutions to true to print U whenever it is used in an assembly() call

Definition at line 396 of file diff_system.h.

Referenced by assembly().

std::vector<Number> 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 1020 of file system.h.

Referenced by ImplicitSystem::adjoint_solve(), SensitivityData::allocate_data(), SensitivityData::allocate_hessian_data(), assemble_qoi(), ExplicitSystem::assemble_qoi(), assemble_qoi_derivative(), ExplicitSystem::assemble_qoi_derivative(), DiffContext::DiffContext(), AdjointResidualErrorEstimator::estimate_error(), FEMContext::reinit(), QoISet::size(), and ImplicitSystem::weighted_sensitivity_adjoint_solve().

NumericVector<Number>* ExplicitSystem::rhs [inherited]

The system matrix. Implicit systems are characterized by the need to solve the linear system Ax=b. This is the right-hand-side vector b.

Definition at line 121 of file explicit_system.h.

Referenced by __libmesh_petsc_diff_solver_residual(), ExplicitSystem::add_system_rhs(), ImplicitSystem::assemble(), ImplicitSystem::assemble_residual_derivatives(), LinearImplicitSystem::assembly(), assembly(), ExplicitSystem::clear(), ContinuationSystem::continuation_solve(), NewtonSolver::line_search(), PetscDiffSolver::solve(), NonlinearImplicitSystem::solve(), NewtonSolver::solve(), LinearImplicitSystem::solve(), FrequencySystem::solve(), ContinuationSystem::solve_tangent(), NewmarkSystem::update_rhs(), ImplicitSystem::weighted_sensitivity_adjoint_solve(), and ImplicitSystem::weighted_sensitivity_solve().

AutoPtr<NumericVector<Number> > System::solution [inherited]

Data structure to hold solution values.

Definition at line 1000 of file system.h.

Referenced by __libmesh_petsc_diff_solver_jacobian(), __libmesh_petsc_diff_solver_residual(), ExactSolution::_compute_error(), UnsteadySolver::advance_timestep(), AdaptiveTimeSolver::advance_timestep(), ContinuationSystem::apply_predictor(), LinearImplicitSystem::assembly(), assembly(), System::clear(), System::compare(), Problem_Interface::computeF(), ContinuationSystem::continuation_solve(), GMVIO::copy_nodal_solution(), UnsteadySolver::du(), DofMap::enforce_constraints_exactly(), PatchRecoveryErrorEstimator::estimate_error(), JumpErrorEstimator::estimate_error(), ExactErrorEstimator::estimate_error(), EigenSystem::get_eigenpair(), System::init_data(), ContinuationSystem::initialize_tangent(), DofMap::max_constraint_error(), mesh_position_get(), ErrorVector::plot_error(), System::project_solution(), System::re_update(), System::read_legacy_data(), System::read_parallel_data(), System::read_serialized_data(), System::reinit(), System::restrict_vectors(), ContinuationSystem::save_current_solution(), VTKIO::solution_to_vtk(), TwostepTimeSolver::solve(), PetscDiffSolver::solve(), NonlinearImplicitSystem::solve(), NewtonSolver::solve(), LinearImplicitSystem::solve(), FrequencySystem::solve(), ContinuationSystem::solve_tangent(), System::update(), System::update_global_solution(), ContinuationSystem::update_solution(), NewmarkSystem::update_u_v_a(), System::write_parallel_data(), and System::write_serialized_data().

Real DifferentiableSystem::time [inherited]

For time-dependent problems, this is the time t at the beginning of the current timestep. But do *not* access this time during an assembly! Use the DiffContext::time value instead to get correct results.

Definition at line 378 of file diff_system.h.

Referenced by UnsteadySolver::advance_timestep(), AdaptiveTimeSolver::advance_timestep(), and TwostepTimeSolver::solve().

AutoPtr<TimeSolver> DifferentiableSystem::time_solver [inherited]

A pointer to the solver object we're going to use. This must be instantiated by the user before solving!

Definition at line 370 of file diff_system.h.

Referenced by assembly(), ContinuationSystem::continuation_solve(), eulerian_residual(), DifferentiableSystem::get_linear_solve_parameters(), DifferentiableSystem::get_linear_solver(), DifferentiableSystem::init_data(), DifferentiableSystem::reinit(), DifferentiableSystem::solve(), and ContinuationSystem::solve_tangent().

bool DifferentiableSystem::use_fixed_solution [inherited]

A boolean to be set to true by systems using elem_fixed_solution, so that the library code will create it. False by default. Warning: if this variable is set to true, it must be before init_data() is called.

Definition at line 430 of file diff_system.h.

Referenced by DifferentiableSystem::clear(), DiffContext::DiffContext(), SteadySolver::element_residual(), EulerSolver::element_residual(), Euler2Solver::element_residual(), FEMContext::reinit(), SteadySolver::side_residual(), EulerSolver::side_residual(), and Euler2Solver::side_residual().

Real FEMSystem::verify_analytic_jacobians

If verify_analytic_jacobian is equal to zero (as it is by default), no numeric jacobians will be calculated unless an overloaded element_time_derivative(), element_constraint(), side_time_derivative(), or side_constraint() function cannot provide an analytic jacobian upon request.

If verify_analytic_jacobian is equal to the positive value tol, then any time a full analytic element jacobian can be calculated it will be tested against a numerical jacobian on the same element, and the program will abort if the relative error (in matrix l1 norms) exceeds tol.

Definition at line 271 of file fem_system.h.

Referenced by assembly().

The documentation for this class was generated from the following files:


Public Types
typedef FEMSystem	sys_type
typedef DifferentiableSystem	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
	FEMSystem (EquationSystems &es, const std::string &name, const unsigned int number)
virtual	~FEMSystem ()
virtual void	clear ()
virtual void	assembly (bool get_residual, bool get_jacobian)
virtual void	solve ()
virtual void	time_evolving (unsigned int var)
virtual void	mesh_x_position (unsigned int sysnum, unsigned int var)
virtual void	mesh_y_position (unsigned int sysnum, unsigned int var)
virtual void	mesh_z_position (unsigned int sysnum, unsigned int var)
void	mesh_position_get ()
void	mesh_position_set ()
virtual bool	eulerian_residual (bool request_jacobian, DiffContext &context)
virtual bool	mass_residual (bool request_jacobian, DiffContext &context)
virtual AutoPtr< DiffContext >	build_context ()
virtual void	init_context (DiffContext &)
virtual void	postprocess ()
virtual void	assemble_qoi (const QoISet &indices=QoISet())
virtual void	assemble_qoi_derivative (const QoISet &indices=QoISet())
virtual void	reinit ()
virtual void	assemble ()
virtual LinearSolver< Number > *	get_linear_solver () const
virtual std::pair< unsigned int, Real >	get_linear_solve_parameters () const
virtual void	release_linear_solver (LinearSolver< Number > *) const
virtual bool	element_time_derivative (bool request_jacobian, DiffContext &)
virtual bool	element_constraint (bool request_jacobian, DiffContext &)
virtual bool	side_time_derivative (bool request_jacobian, DiffContext &)
virtual bool	side_constraint (bool request_jacobian, DiffContext &)
virtual void	element_postprocess (DiffContext &)
virtual void	side_postprocess (DiffContext &)
virtual void	element_qoi (DiffContext &, const QoISet &)
virtual void	element_qoi_derivative (DiffContext &, const QoISet &)
virtual void	side_qoi (DiffContext &, const QoISet &)
virtual void	side_qoi_derivative (DiffContext &, const QoISet &)
virtual bool	side_mass_residual (bool request_jacobian, DiffContext &)
sys_type &	system ()
virtual std::string	system_type () 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)
unsigned int	n_matrices () const
void	init ()
virtual void	update ()
virtual void	qoi_parameter_sensitivity (const QoISet &qoi_indices, const ParameterVector &parameters, SensitivityData &sensitivities)
virtual bool	compare (const System &other_system, const Real threshold, const bool verbose) const
const std::string &	name () const
void	project_solution (Number fptr(const Point &p, const Parameters &parameters, const std::string &sys_name, const std::string &unknown_name), Gradient gptr(const Point &p, const Parameters &parameters, const std::string &sys_name, const std::string &unknown_name), Parameters &parameters) const
void	project_vector (Number fptr(const Point &p, const Parameters &parameters, const std::string &sys_name, const std::string &unknown_name), Gradient gptr(const Point &p, const Parameters &parameters, const std::string &sys_name, const std::string &unknown_name), Parameters &parameters, NumericVector< Number > &new_vector) const
unsigned int	number () const
void	update_global_solution (std::vector< Number > &global_soln) const
void	update_global_solution (std::vector< Number > &global_soln, const unsigned int dest_proc) const
const MeshBase &	get_mesh () const
MeshBase &	get_mesh ()
const DofMap &	get_dof_map () const
DofMap &	get_dof_map ()
const EquationSystems &	get_equation_systems () const
EquationSystems &	get_equation_systems ()
bool	active () const
void	activate ()
void	deactivate ()
vectors_iterator	vectors_begin ()
const_vectors_iterator	vectors_begin () const
vectors_iterator	vectors_end ()
const_vectors_iterator	vectors_end () const
NumericVector< Number > &	add_vector (const std::string &vec_name, const bool projections=true)
bool &	project_solution_on_reinit (void)
bool	have_vector (const std::string &vec_name) const
const NumericVector< Number > *	request_vector (const std::string &vec_name) const
NumericVector< Number > *	request_vector (const std::string &vec_name)
const NumericVector< Number > *	request_vector (const unsigned int vec_num) const
NumericVector< Number > *	request_vector (const unsigned int vec_num)
const NumericVector< Number > &	get_vector (const std::string &vec_name) const
NumericVector< Number > &	get_vector (const std::string &vec_name)
const NumericVector< Number > &	get_vector (const unsigned int vec_num) const
NumericVector< Number > &	get_vector (const unsigned int vec_num)
const std::string &	vector_name (const unsigned int vec_num)
NumericVector< Number > &	add_adjoint_solution (unsigned int i=0)
NumericVector< Number > &	get_adjoint_solution (unsigned int i=0)
const NumericVector< Number > &	get_adjoint_solution (unsigned int i=0) const
NumericVector< Number > &	add_sensitivity_solution (unsigned int i=0)
NumericVector< Number > &	get_sensitivity_solution (unsigned int i=0)
const NumericVector< Number > &	get_sensitivity_solution (unsigned int i=0) const
NumericVector< Number > &	add_weighted_sensitivity_adjoint_solution (unsigned int i=0)
NumericVector< Number > &	get_weighted_sensitivity_adjoint_solution (unsigned int i=0)
const NumericVector< Number > &	get_weighted_sensitivity_adjoint_solution (unsigned int i=0) const
NumericVector< Number > &	add_weighted_sensitivity_solution ()
NumericVector< Number > &	get_weighted_sensitivity_solution ()
const NumericVector< Number > &	get_weighted_sensitivity_solution () const
NumericVector< Number > &	add_adjoint_rhs (unsigned int i=0)
NumericVector< Number > &	get_adjoint_rhs (unsigned int i=0)
const NumericVector< Number > &	get_adjoint_rhs (unsigned int i=0) const
NumericVector< Number > &	add_sensitivity_rhs (unsigned int i=0)
NumericVector< Number > &	get_sensitivity_rhs (unsigned int i=0)
const NumericVector< Number > &	get_sensitivity_rhs (unsigned int i=0) const
unsigned int	n_vectors () const
unsigned int	n_vars () const
unsigned int	n_dofs () const
unsigned int	n_active_dofs () const
unsigned int	n_constrained_dofs () const
unsigned int	n_local_dofs () const
unsigned int	add_variable (const std::string &var, const FEType &type, const std::set< subdomain_id_type > *const active_subdomains=NULL)
unsigned int	add_variable (const std::string &var, const Order order=FIRST, const FEFamily=LAGRANGE, const std::set< subdomain_id_type > *const active_subdomains=NULL)
const Variable &	variable (unsigned int var) const
bool	has_variable (const std::string &var) const
const std::string &	variable_name (const unsigned int i) const
unsigned short int	variable_number (const std::string &var) const
const FEType &	variable_type (const unsigned int i) const
const FEType &	variable_type (const std::string &var) const
Real	calculate_norm (NumericVector< Number > &v, unsigned int var=0, FEMNormType norm_type=L2) const
Real	calculate_norm (NumericVector< Number > &v, const SystemNorm &norm) const
void	read_header (Xdr &io, const std::string &version, const bool read_header=true, const bool read_additional_data=true, const bool read_legacy_format=false)
void	read_legacy_data (Xdr &io, const bool read_additional_data=true)
void	read_serialized_data (Xdr &io, const bool read_additional_data=true)
void	read_parallel_data (Xdr &io, const bool read_additional_data)
void	write_header (Xdr &io, const std::string &version, const bool write_additional_data) const
void	write_serialized_data (Xdr &io, const bool write_additional_data=true) const
void	write_parallel_data (Xdr &io, const bool write_additional_data) const
std::string	get_info () const
void	attach_init_function (void fptr(EquationSystems &es, const std::string &name))
void	attach_assemble_function (void fptr(EquationSystems &es, const std::string &name))
void	attach_constraint_function (void fptr(EquationSystems &es, const std::string &name))
void	attach_QOI_function (void fptr(EquationSystems &es, const std::string &name, const QoISet &qoi_indices))
void	attach_QOI_derivative (void fptr(EquationSystems &es, const std::string &name, const QoISet &qoi_indices))
virtual void	user_initialization ()
virtual void	user_assembly ()
virtual void	user_constrain ()
virtual void	user_QOI (const QoISet &qoi_indices)
virtual void	user_QOI_derivative (const QoISet &qoi_indices)
virtual void	re_update ()
virtual void	restrict_vectors ()
virtual void	prolong_vectors ()
Number	current_solution (const unsigned int global_dof_number) const
void	local_dof_indices (const unsigned int var, std::set< unsigned int > &var_indices) const
void	zero_variable (NumericVector< Number > &v, unsigned int var_num) const
Static Public Member Functions
static std::string	get_info ()
static void	print_info ()
static unsigned int	n_objects ()
Public Attributes
bool	fe_reinit_during_postprocess
int	extra_quadrature_order
Real	numerical_jacobian_h
Real	verify_analytic_jacobians
bool	compute_internal_sides
bool	postprocess_sides
bool	assemble_qoi_sides
AutoPtr< TimeSolver >	time_solver
Real	time
Real	deltat
bool	print_solution_norms
bool	print_solutions
bool	print_residual_norms
bool	print_residuals
bool	print_jacobian_norms
bool	print_jacobians
bool	print_element_jacobians
bool	use_fixed_solution
SparseMatrix< Number > *	matrix
NumericVector< Number > *	rhs
bool	assemble_before_solve
AutoPtr< NumericVector< Number > >	solution
AutoPtr< NumericVector< Number > >	current_local_solution
std::vector< Number >	qoi
Protected Types
typedef bool(TimeSolver::*	TimeSolverResPtr )(bool, DiffContext &)
typedef std::map< std::string, std::pair< unsigned int, unsigned int > >	Counts
Protected Member Functions
virtual void	init_data ()
void	numerical_jacobian (TimeSolverResPtr res, FEMContext &context)
void	numerical_elem_jacobian (FEMContext &context)
void	numerical_side_jacobian (FEMContext &context)
virtual void	init_matrices ()
void	project_vector (NumericVector< Number > &) const
void	project_vector (const NumericVector< Number > &, NumericVector< Number > &) const
void	increment_constructor_count (const std::string &name)
void	increment_destructor_count (const std::string &name)
Protected Attributes
unsigned int	_mesh_sys
unsigned int	_mesh_x_var
unsigned int	_mesh_y_var
unsigned int	_mesh_z_var
std::vector< bool >	_time_evolving
Static Protected Attributes
static Counts	_counts
static Threads::atomic < unsigned int >	_n_objects
static Threads::spin_mutex	_mutex