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_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 256 of file implicit_system.h.

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

Definition at line 405 of file system.h.

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

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

Definition at line 105 of file reference_counter.h.

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

Matrix iterator typedefs.

Definition at line 255 of file implicit_system.h.

typedef 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 404 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 1374 of file system.h.

References System::_active.

01375 {
01376   _active = true;
01377 }

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

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

Definition at line 1366 of file system.h.

References System::_active.

01367 {
01368   return _active;  
01369 }

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 1382 of file system.h.

References System::_active.

01383 {
01384   _active = false;
01385 }

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 1358 of file system.h.

References System::_dof_map.

01359 {
01360   return *_dof_map;
01361 }

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

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

Definition at line 1350 of file system.h.

References System::_dof_map.

Referenced by __libmesh_petsc_diff_solver_jacobian(), __libmesh_petsc_diff_solver_residual(), ExactSolution::_compute_error(), HPCoarsenTest::add_projection(), ImplicitSystem::adjoint_solve(), 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().

01351 {
01352   return *_dof_map;
01353 }

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

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

Definition at line 383 of file system.h.

References System::_equation_systems.

00383 { return _equation_systems; }

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

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

Definition at line 378 of file system.h.

References System::_equation_systems.

00378 { 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 1342 of file system.h.

References System::_mesh.

01343 {
01344   return _mesh;
01345 }

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

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

Definition at line 1334 of file system.h.

References System::_mesh.

Referenced by ExactSolution::_compute_error(), HPCoarsenTest::add_projection(), 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().

01335 {
01336   return _mesh;
01337 }

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 357 of file implicit_system.h.

References ImplicitSystem::_matrices.

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

00358 {
00359   return (_matrices.count(mat_name));
00360 }

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

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

Definition at line 1444 of file system.h.

References System::_vectors.

Referenced by System::add_vector().

01445 {
01446   return (_vectors.count(vec_name));
01447 }

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 1436 of file system.h.

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

01437 {
01438   return this->n_dofs() - this->n_constrained_dofs();
01439 }

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 364 of file implicit_system.h.

References ImplicitSystem::_matrices.

00365 {
00366  return _matrices.size(); 
00367 }

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 1390 of file system.h.

References System::_variables.

Referenced by System::add_variable(), EquationSystems::build_discontinuous_solution_vector(), EquationSystems::build_solution_vector(), System::calculate_norm(), DiffContext::DiffContext(), JumpErrorEstimator::estimate_error(), ExactErrorEstimator::estimate_error(), ErrorEstimator::estimate_errors(), 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().

01391 {
01392   return _variables.size();
01393 }

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

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

Definition at line 1452 of file system.h.

References System::_vectors.

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

01453 {
01454   return _vectors.size(); 
01455 }

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

Returns:: the system name.

Definition at line 1318 of file system.h.

References System::_sys_name.

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

01319 {
01320   return _sys_name;
01321 }

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

Returns:: the system number.

Definition at line 1326 of file system.h.

References System::_sys_number.

Referenced by System::add_variable(), 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().

01327 {
01328   return _sys_number;
01329 }

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 445 of file system.h.

References System::_solution_projection.

00446     { 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


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