LinearImplicitSystem Class Reference

#include <linear_implicit_system.h>

Inheritance diagram for LinearImplicitSystem:

List of all members.

Public Types

typedef LinearImplicitSystem sys_type

typedef ImplicitSystem Parent

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

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

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

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

Public Member Functions

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

virtual ~LinearImplicitSystem ()

sys_type & system ()

virtual void clear ()

virtual void reinit ()

virtual void assemble ()

virtual void solve ()

virtual LinearSolver< Number > * get_linear_solver () const

virtual void release_linear_solver (LinearSolver< Number > *) const

virtual void assembly (bool get_residual, bool get_jacobian)

virtual std::string system_type () const

unsigned int n_linear_iterations () const

Real final_linear_residual () const

void attach_shell_matrix (ShellMatrix< Number > *shell_matrix)

void detach_shell_matrix (void)

ShellMatrix< Number > * get_shell_matrix (void)

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

virtual void assemble_residual_derivatives (const ParameterVector &parameters)

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

virtual std::pair< unsigned
int, Real > weighted_sensitivity_solve (const ParameterVector &parameters, const ParameterVector &weights)

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

virtual std::pair< unsigned
int, Real > weighted_sensitivity_adjoint_solve (const ParameterVector &parameters, const ParameterVector &weights, const QoISet &qoi_indices=QoISet())

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

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

virtual void qoi_parameter_hessian (const QoISet &qoi_indices, const ParameterVector &parameters, SensitivityData &hessian)

virtual void qoi_parameter_hessian_vector_product (const QoISet &qoi_indices, const ParameterVector &parameters, const ParameterVector &vector, SensitivityData &product)

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

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

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

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

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

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

unsigned int n_matrices () const

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

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

void init ()

virtual void update ()

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

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

const std::string & name () const

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

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

unsigned int number () const

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

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

const MeshBase & get_mesh () const

MeshBase & get_mesh ()

const DofMap & get_dof_map () const

DofMap & get_dof_map ()

const EquationSystems & get_equation_systems () const

EquationSystems & get_equation_systems ()

bool active () const

void activate ()

void deactivate ()

vectors_iterator vectors_begin ()

const_vectors_iterator vectors_begin () const

vectors_iterator vectors_end ()

const_vectors_iterator vectors_end () const

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

bool & project_solution_on_reinit (void)

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

NumericVector< Number > & add_weighted_sensitivity_solution ()

NumericVector< Number > & get_weighted_sensitivity_solution ()

const NumericVector< Number > & get_weighted_sensitivity_solution () const

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

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

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

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

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

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

unsigned int n_vectors () const

unsigned int n_vars () const

unsigned int n_dofs () const

unsigned int n_active_dofs () const

unsigned int n_constrained_dofs () const

unsigned int n_local_dofs () const

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

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

const Variable & variable (unsigned int var) const

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

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

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

const FEType & variable_type (const unsigned int i) const

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

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

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

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

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

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

void read_parallel_data (Xdr &io, const bool read_additional_data)

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

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

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

std::string get_info () const

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

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

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

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

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

virtual void user_initialization ()

virtual void user_assembly ()

virtual void user_constrain ()

virtual void user_QOI (const QoISet &qoi_indices)

virtual void user_QOI_derivative (const QoISet &qoi_indices)

virtual void re_update ()

virtual void restrict_vectors ()

virtual void prolong_vectors ()

Number current_solution (const unsigned int global_dof_number) const

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

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

Static Public Member Functions

static std::string get_info ()

static void print_info ()

static unsigned int n_objects ()

Public Attributes

AutoPtr< LinearSolver< Number > > linear_solver

SparseMatrix< Number > * matrix

NumericVector< Number > * rhs

bool assemble_before_solve

AutoPtr< NumericVector< Number > > solution

AutoPtr< NumericVector< Number > > current_local_solution

std::vector< Number > qoi

Protected Types

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

Protected Member Functions

virtual void init_data ()

virtual void init_matrices ()

void project_vector (NumericVector< Number > &) const

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

void increment_constructor_count (const std::string &name)

void increment_destructor_count (const std::string &name)

Protected Attributes

unsigned int _n_linear_iterations

Real _final_linear_residual

ShellMatrix< Number > * _shell_matrix

Static Protected Attributes

static Counts _counts

static Threads::atomic
< unsigned int > _n_objects

static Threads::spin_mutex _mutex

Detailed Description

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

Definition at line 47 of file linear_implicit_system.h.

Member Typedef Documentation

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

Definition at line 268 of file implicit_system.h.

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

Definition at line 411 of file system.h.

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

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

Definition at line 105 of file reference_counter.h.

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

Matrix iterator typedefs.

Definition at line 267 of file implicit_system.h.

typedef ImplicitSystem LinearImplicitSystem::Parent

The type of the parent.

Reimplemented from ImplicitSystem.

Definition at line 72 of file linear_implicit_system.h.

typedef LinearImplicitSystem LinearImplicitSystem::sys_type

The type of system.

Reimplemented from ImplicitSystem.

Reimplemented in NewmarkSystem.

Definition at line 67 of file linear_implicit_system.h.

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

Vector iterator typedefs.

Definition at line 410 of file system.h.

Constructor & Destructor Documentation

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

Constructor. Optionally initializes required data structures.

Definition at line 35 of file linear_implicit_system.C.

00037                                                                        :
00038   
00039   Parent                 (es, name, number),
00040   linear_solver          (LinearSolver<Number>::build()),
00041   _n_linear_iterations   (0),
00042   _final_linear_residual (1.e20),
00043   _shell_matrix(NULL)
00044 {
00045 }

LinearImplicitSystem::~LinearImplicitSystem ( ) [virtual]

Destructor.

Definition at line 49 of file linear_implicit_system.C.

References clear().

00050 {
00051   // Clear data
00052   this->clear();
00053 }

Member Function Documentation

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

Activates the system. Only active systems are solved.

Definition at line 1380 of file system.h.

References System::_active.

01381 {
01382   _active = true;
01383 }

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

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

Definition at line 1372 of file system.h.

References System::_active.

01373 {
01374   return _active;  
01375 }

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

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

Definition at line 824 of file system.C.

References System::add_vector().

Referenced by FEMSystem::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 }

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

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

Reimplemented from ImplicitSystem.

Reimplemented in FrequencySystem, and NewmarkSystem.

Definition at line 97 of file linear_implicit_system.h.

Referenced by assembly(), and solve().

00097 { ImplicitSystem::assemble(); }

void ExplicitSystem::assemble_qoi ( const QoISet & qoi_indices = QoISet() ) [virtual, inherited]

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

Reimplemented from System.

Reimplemented in FEMSystem.

Definition at line 89 of file explicit_system.C.

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

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

void ExplicitSystem::assemble_qoi_derivative ( const QoISet & qoi_indices = QoISet() ) [virtual, inherited]

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

Reimplemented from System.

Reimplemented in FEMSystem.

Definition at line 102 of file explicit_system.C.

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

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

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

void 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 LinearImplicitSystem::assembly	(	bool	get_residual,
		bool	get_jacobian
	)			`[virtual]`

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

Reimplemented from ImplicitSystem.

Definition at line 333 of file linear_implicit_system.C.

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

00335 {
00336   // Residual R(u(p),p) := A(p)*u(p) - b(p)
00337   // partial R / partial u = A
00338 
00339   this->assemble();
00340   this->rhs->close();
00341   this->matrix->close();
00342 
00343   *(this->rhs) *= -1.0;
00344   this->rhs->add_vector(*this->solution, *this->matrix);
00345 }

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

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

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

Definition at line 125 of file linear_implicit_system.C.

References _shell_matrix.

Referenced by detach_shell_matrix().

00126 {
00127   _shell_matrix = shell_matrix;
00128 }

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 LinearImplicitSystem::clear ( ) [virtual]

Clear all the data structures associated with the system.

Reimplemented from ImplicitSystem.

Reimplemented in FrequencySystem, and NewmarkSystem.

Definition at line 57 of file linear_implicit_system.C.

References ImplicitSystem::clear(), and linear_solver.

Referenced by ~LinearImplicitSystem().

00058 {
00059   // clear the linear solver
00060   linear_solver->clear();
00061   
00062   // clear the parent data
00063   Parent::clear();
00064 }

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(), FEMSystem::eulerian_residual(), PatchRecoveryErrorEstimator::EstimateError::operator()(), FEMContext::reinit(), HPCoarsenTest::select_refinement(), VTKIO::solution_to_vtk(), EnsightIO::write_scalar_ascii(), and EnsightIO::write_vector_ascii().

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

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

Deactivates the system. Only active systems are solved.

Definition at line 1388 of file system.h.

References System::_active.

01389 {
01390   _active = false;
01391 }

void LinearImplicitSystem::detach_shell_matrix ( void ) [inline]

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

Definition at line 161 of file linear_implicit_system.h.

References attach_shell_matrix().

00161 { attach_shell_matrix(NULL); }

Real LinearImplicitSystem::final_linear_residual ( ) const [inline]

Returns the final residual for the linear system solve.

Definition at line 145 of file linear_implicit_system.h.

References _final_linear_residual.

00145 { return _final_linear_residual; }

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

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

Uses the forward sensitivity method.

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

Reimplemented from System.

Definition at line 749 of file implicit_system.C.

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

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(), FEMSystem::assemble_qoi_derivative(), and ImplicitSystem::weighted_sensitivity_adjoint_solve().

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

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

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

Definition at line 784 of file system.C.

References System::get_vector().

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

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

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

Definition at line 774 of file system.C.

References System::get_vector().

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

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

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

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

Definition at line 1364 of file system.h.

References System::_dof_map.

01365 {
01366   return *_dof_map;
01367 }

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

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

Definition at line 1356 of file system.h.

References System::_dof_map.

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

01357 {
01358   return *_dof_map;
01359 }

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

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

Definition at line 389 of file system.h.

References System::_equation_systems.

00389 { return _equation_systems; }

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

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

Definition at line 384 of file system.h.

References System::_equation_systems.

00384 { return _equation_systems; }

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

Gets a string containing the reference information.

Definition at line 45 of file reference_counter.C.

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

Referenced by ReferenceCounter::print_info().

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

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

Returns:: a string containing information about the system.

Definition at line 1233 of file system.C.

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

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

std::pair< unsigned int, Real > ImplicitSystem::get_linear_solve_parameters ( ) const [virtual, inherited]

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

Reimplemented in DifferentiableSystem, and NonlinearImplicitSystem.

Definition at line 1261 of file implicit_system.C.

References System::get_equation_systems().

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

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

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

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

Reimplemented from ImplicitSystem.

Definition at line 320 of file linear_implicit_system.C.

References linear_solver.

00321 {
00322   return linear_solver.get();
00323 }

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

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

Definition at line 245 of file implicit_system.C.

References ImplicitSystem::_matrices.

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

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

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

Definition at line 226 of file implicit_system.C.

References ImplicitSystem::_matrices.

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

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

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

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

Definition at line 1348 of file system.h.

References System::_mesh.

01349 {
01350   return _mesh;
01351 }

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

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

Definition at line 1340 of file system.h.

References System::_mesh.

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

01341 {
01342   return _mesh;
01343 }

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

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

Definition at line 874 of file system.C.

References System::get_vector().

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

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

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

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

Definition at line 864 of file system.C.

References System::get_vector().

Referenced by ImplicitSystem::sensitivity_solve().

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

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

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

Definition at line 733 of file system.C.

References System::get_vector().

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

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

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

Definition at line 723 of file system.C.

References System::get_vector().

Referenced by ImplicitSystem::sensitivity_solve().

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

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

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

Definition at line 167 of file linear_implicit_system.h.

References _shell_matrix.

00167 { return _shell_matrix; }

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

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

Definition at line 682 of file system.C.

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

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

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

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

Definition at line 666 of file system.C.

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

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

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

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

Definition at line 647 of file system.C.

References System::_vectors.

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

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

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

Definition at line 628 of file system.C.

References System::_vectors.

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

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

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

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

Definition at line 814 of file system.C.

References System::get_vector().

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

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

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

Definition at line 804 of file system.C.

References System::get_vector().

Referenced by ImplicitSystem::weighted_sensitivity_adjoint_solve().

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

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

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

Definition at line 757 of file system.C.

References System::get_vector().

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

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

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

Definition at line 750 of file system.C.

References System::get_vector().

Referenced by ImplicitSystem::weighted_sensitivity_solve().

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

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

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

Definition at line 935 of file system.C.

References System::_variable_numbers.

Referenced by GMVIO::copy_nodal_solution().

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

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

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

Definition at line 369 of file implicit_system.h.

References ImplicitSystem::_matrices.

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

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

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

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

Definition at line 1450 of file system.h.

References System::_vectors.

Referenced by System::add_vector().

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

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

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

Definition at line 149 of file reference_counter.h.

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

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

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

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

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

Definition at line 167 of file reference_counter.h.

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

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

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

void System::init ( ) [inherited]

Initializes degrees of freedom on the current mesh. Sets the

Definition at line 158 of file system.C.

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

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

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

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

Reimplemented from ExplicitSystem.

Reimplemented in ContinuationSystem, DifferentiableSystem, FEMSystem, and FrequencySystem.

Definition at line 85 of file implicit_system.C.

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

Referenced by DifferentiableSystem::init_data().

00086 {
00087   // initialize parent data
00088   Parent::init_data();
00089 
00090   // Add the system matrix.
00091   this->add_system_matrix ();
00092   
00093   // Initialize the matrices for the system
00094   this->init_matrices ();
00095 }

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 }

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

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

Definition at line 1442 of file system.h.

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

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

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

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

Definition at line 95 of file system.C.

References System::_dof_map.

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

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

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

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

Definition at line 88 of file system.C.

References System::_dof_map.

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

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

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

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

Definition at line 140 of file linear_implicit_system.h.

References _n_linear_iterations.

00140 { return _n_linear_iterations; }

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

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

Definition at line 110 of file system.C.

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

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

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

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

Returns:: the number of matrices handled by this system

Definition at line 376 of file implicit_system.h.

References ImplicitSystem::_matrices.

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

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

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

Definition at line 76 of file reference_counter.h.

References ReferenceCounter::_n_objects.

00077   { return _n_objects; }

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

Returns:: the number of variables in the system

Definition at line 1396 of file system.h.

References System::_variables.

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

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

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

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

Definition at line 1458 of file system.h.

References System::_vectors.

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

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

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

Returns:: the system name.

Definition at line 1324 of file system.h.

References System::_sys_name.

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

01325 {
01326   return _sys_name;
01327 }

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

Returns:: the system number.

Definition at line 1332 of file system.h.

References System::_sys_number.

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

01333 {
01334   return _sys_number;
01335 }

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

Prints the reference information to std::cout.

Definition at line 83 of file reference_counter.C.

References ReferenceCounter::get_info().

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

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

Projects the continuous functions onto the current solution.

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

Definition at line 237 of file system_projection.C.

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

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

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

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

Definition at line 451 of file system.h.

References System::_solution_projection.

00452     { return _solution_projection; }

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

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

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

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

Definition at line 60 of file system_projection.C.

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

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

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

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

Definition at line 43 of file system_projection.C.

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

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

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

Projects the continuous functions onto the current mesh.

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

Definition at line 258 of file system_projection.C.

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

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

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

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

Prolong vectors after the mesh has refined

Definition at line 255 of file system.C.

References System::restrict_vectors().

Referenced by EquationSystems::reinit().

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

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

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

Reimplemented from System.

Definition at line 1004 of file implicit_system.C.

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

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

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

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

Reimplemented from System.

Definition at line 825 of file implicit_system.C.

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

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

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

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

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

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

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

Definition at line 383 of file system.C.

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

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

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

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

Definition at line 314 of file system.C.

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

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

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

Reads the basic data header for this System.

Definition at line 78 of file system_io.C.

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

Referenced by EquationSystems::_read_impl().

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

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

Reads additional data, namely vectors, for this System.

Definition at line 241 of file system_io.C.

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

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

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

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

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

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

for each additional vector in the object

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

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

Definition at line 419 of file system_io.C.

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

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

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

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

Definition at line 559 of file system_io.C.

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

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

void LinearImplicitSystem::reinit ( ) [virtual]

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

Reimplemented from ImplicitSystem.

Reimplemented in NewmarkSystem.

Definition at line 68 of file linear_implicit_system.C.

References linear_solver, and ImplicitSystem::reinit().

00069 {
00070   // re-initialize the linear solver interface
00071   linear_solver->clear();
00072   
00073   // initialize parent data
00074   Parent::reinit();  
00075 }

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

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

Reimplemented from ImplicitSystem.

Definition at line 327 of file linear_implicit_system.C.

00328 {
00329 }

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

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

Definition at line 213 of file implicit_system.C.

References ImplicitSystem::_matrices.

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

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

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

Definition at line 200 of file implicit_system.C.

References ImplicitSystem::_matrices.

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

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

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

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

Definition at line 611 of file system.C.

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

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

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

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

Definition at line 594 of file system.C.

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

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

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

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

Definition at line 582 of file system.C.

References System::_vectors.

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

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

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

Definition at line 570 of file system.C.

References System::_vectors.

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

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

Restrict vectors after the mesh has coarsened

Definition at line 219 of file system.C.

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

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

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

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

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

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

Reimplemented from System.

Definition at line 282 of file implicit_system.C.

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

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

void LinearImplicitSystem::solve ( ) [virtual]

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

Reimplemented from ExplicitSystem.

Reimplemented in FrequencySystem.

Definition at line 79 of file linear_implicit_system.C.

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

00080 {
00081   if (this->assemble_before_solve)
00082     // Assemble the linear system
00083     this->assemble (); 
00084 
00085   // Log how long the linear solve takes.
00086   // This gets done by the LinearSolver classes now [RHS]
00087   // START_LOG("solve()", "System");
00088 
00089   // Get a reference to the EquationSystems
00090   const EquationSystems& es =
00091     this->get_equation_systems();
00092   
00093   // Get the user-specifiied linear solver tolerance
00094   const Real tol            =
00095     es.parameters.get<Real>("linear solver tolerance");
00096 
00097   // Get the user-specified maximum # of linear solver iterations
00098   const unsigned int maxits =
00099     es.parameters.get<unsigned int>("linear solver maximum iterations");
00100 
00101   // Solve the linear system.  Several cases:
00102   std::pair<unsigned int, Real> rval = std::make_pair(0,0.0);
00103   if(_shell_matrix)
00104     // 1.) Shell matrix with or without user-supplied preconditioner.
00105     rval = linear_solver->solve(*_shell_matrix, this->request_matrix("Preconditioner"), *solution, *rhs, tol, maxits);
00106   else
00107     // 2.) No shell matrix, with or without user-supplied preconditioner
00108     rval = linear_solver->solve (*matrix, this->request_matrix("Preconditioner"), *solution, *rhs, tol, maxits);
00109 
00110   // Store the number of linear iterations required to
00111   // solve and the final residual.
00112   _n_linear_iterations   = rval.first;
00113   _final_linear_residual = rval.second;
00114     
00115   // Stop logging the linear solve
00116   // This gets done by the LinearSolver classes now [RHS]
00117   // STOP_LOG("solve()", "System");
00118 
00119   // Update the system after the solve
00120   this->update();  
00121 }

sys_type& LinearImplicitSystem::system ( ) [inline]

Returns:: a clever pointer to the system.

Reimplemented from ImplicitSystem.

Definition at line 77 of file linear_implicit_system.h.

00077 { return *this; }

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

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

Reimplemented from ImplicitSystem.

Reimplemented in FrequencySystem, and NewmarkSystem.

Definition at line 126 of file linear_implicit_system.h.

00126 { return "LinearImplicit"; }

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

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

Definition at line 295 of file system.C.

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

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

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

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

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

Definition at line 539 of file system.C.

References System::solution.

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

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

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

Definition at line 530 of file system.C.

References System::solution.

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

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

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

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

Definition at line 1367 of file system.C.

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

Referenced by System::assemble().

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

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

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

Definition at line 1377 of file system.C.

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

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

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

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

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

Definition at line 1358 of file system.C.

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

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

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

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

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

Definition at line 1387 of file system.C.

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

Referenced by System::assemble_qoi().

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

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

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

Definition at line 1397 of file system.C.

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

Referenced by System::assemble_qoi_derivative().

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

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

Return a constant reference to Variable var.

Definition at line 1404 of file system.h.

References System::_variables.

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

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

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

Returns:: the name of variable i.

Definition at line 1414 of file system.h.

References System::_variables.

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

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

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

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

Definition at line 942 of file system.C.

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

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

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

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

Returns:: the finite element type for variable var.

Definition at line 1434 of file system.h.

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

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

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

Returns:: the finite element type variable number i.

Definition at line 1424 of file system.h.

References System::_variables.

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

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

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

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

Definition at line 698 of file system.C.

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

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

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

Beginning of vectors container

Definition at line 1470 of file system.h.

References System::_vectors.

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

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

Beginning of vectors container

Definition at line 1464 of file system.h.

References System::_vectors.

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

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

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

End of vectors container

Definition at line 1482 of file system.h.

References System::_vectors.

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

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

End of vectors container

Definition at line 1476 of file system.h.

References System::_vectors.

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

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

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

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

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

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

Reimplemented from System.

Definition at line 404 of file implicit_system.C.

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

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

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

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

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

Reimplemented from System.

Definition at line 540 of file implicit_system.C.

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

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

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

Writes the basic data header for this System.

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

for this system

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

for each variable in the system

6.) The name of the variable (string)

7.) Combined in an FEType:

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

end variable loop

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

for each additional vector in the system object

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

end system

Definition at line 815 of file system_io.C.

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

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

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

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

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

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

for each additional vector in the object

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

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

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

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

for each additional vector in the object

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

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

Definition at line 1213 of file system_io.C.

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

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

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

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

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

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

for each additional vector in the object

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

Definition at line 1370 of file system_io.C.

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

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

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

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

Definition at line 995 of file system.C.

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

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

Member Data Documentation

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

Actually holds the data.

Definition at line 110 of file reference_counter.h.

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

Real LinearImplicitSystem::_final_linear_residual [protected]

The final residual for the linear system Ax=b.

Definition at line 180 of file linear_implicit_system.h.

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

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

Mutual exclusion object to enable thread-safe reference counting.

Definition at line 123 of file reference_counter.h.

unsigned int LinearImplicitSystem::_n_linear_iterations [protected]

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

Definition at line 175 of file linear_implicit_system.h.

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

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

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

Definition at line 118 of file reference_counter.h.

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

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

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

Definition at line 185 of file linear_implicit_system.h.

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

bool System::assemble_before_solve [inherited]

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

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

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

Definition at line 985 of file system.h.

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

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

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

Definition at line 1012 of file system.h.

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

AutoPtr<LinearSolver<Number> > LinearImplicitSystem::linear_solver

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

Definition at line 134 of file linear_implicit_system.h.

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

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

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

Definition at line 326 of file implicit_system.h.

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

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

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

Definition at line 1020 of file system.h.

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

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

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

Definition at line 121 of file explicit_system.h.

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

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

Data structure to hold solution values.

Definition at line 1000 of file system.h.

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

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


Public Types
typedef LinearImplicitSystem	sys_type
typedef ImplicitSystem	Parent
typedef std::map< std::string, SparseMatrix< Number > * >::iterator	matrices_iterator
typedef std::map< std::string, SparseMatrix< Number > * >::const_iterator	const_matrices_iterator
typedef std::map< std::string, NumericVector< Number > * >::iterator	vectors_iterator
typedef std::map< std::string, NumericVector< Number > * >::const_iterator	const_vectors_iterator
Public Member Functions
	LinearImplicitSystem (EquationSystems &es, const std::string &name, const unsigned int number)
virtual	~LinearImplicitSystem ()
sys_type &	system ()
virtual void	clear ()
virtual void	reinit ()
virtual void	assemble ()
virtual void	solve ()
virtual LinearSolver< Number > *	get_linear_solver () const
virtual void	release_linear_solver (LinearSolver< Number > *) const
virtual void	assembly (bool get_residual, bool get_jacobian)
virtual std::string	system_type () const
unsigned int	n_linear_iterations () const
Real	final_linear_residual () const
void	attach_shell_matrix (ShellMatrix< Number > *shell_matrix)
void	detach_shell_matrix (void)
ShellMatrix< Number > *	get_shell_matrix (void)
virtual std::pair< unsigned int, Real >	get_linear_solve_parameters () const
virtual void	assemble_residual_derivatives (const ParameterVector &parameters)
virtual std::pair< unsigned int, Real >	sensitivity_solve (const ParameterVector &parameters)
virtual std::pair< unsigned int, Real >	weighted_sensitivity_solve (const ParameterVector &parameters, const ParameterVector &weights)
virtual std::pair< unsigned int, Real >	adjoint_solve (const QoISet &qoi_indices=QoISet())
virtual std::pair< unsigned int, Real >	weighted_sensitivity_adjoint_solve (const ParameterVector &parameters, const ParameterVector &weights, const QoISet &qoi_indices=QoISet())
virtual void	adjoint_qoi_parameter_sensitivity (const QoISet &qoi_indices, const ParameterVector &parameters, SensitivityData &sensitivities)
virtual void	forward_qoi_parameter_sensitivity (const QoISet &qoi_indices, const ParameterVector &parameters, SensitivityData &sensitivities)
virtual void	qoi_parameter_hessian (const QoISet &qoi_indices, const ParameterVector &parameters, SensitivityData &hessian)
virtual void	qoi_parameter_hessian_vector_product (const QoISet &qoi_indices, const ParameterVector &parameters, const ParameterVector &vector, SensitivityData &product)
SparseMatrix< Number > &	add_matrix (const std::string &mat_name)
bool	have_matrix (const std::string &mat_name) const
const SparseMatrix< Number > *	request_matrix (const std::string &mat_name) const
SparseMatrix< Number > *	request_matrix (const std::string &mat_name)
const SparseMatrix< Number > &	get_matrix (const std::string &mat_name) const
SparseMatrix< Number > &	get_matrix (const std::string &mat_name)
unsigned int	n_matrices () const
virtual void	assemble_qoi (const QoISet &qoi_indices=QoISet())
virtual void	assemble_qoi_derivative (const QoISet &qoi_indices=QoISet())
void	init ()
virtual void	update ()
virtual void	qoi_parameter_sensitivity (const QoISet &qoi_indices, const ParameterVector &parameters, SensitivityData &sensitivities)
virtual bool	compare (const System &other_system, const Real threshold, const bool verbose) const
const std::string &	name () const
void	project_solution (Number fptr(const Point &p, const Parameters &parameters, const std::string &sys_name, const std::string &unknown_name), Gradient gptr(const Point &p, const Parameters &parameters, const std::string &sys_name, const std::string &unknown_name), Parameters &parameters) const
void	project_vector (Number fptr(const Point &p, const Parameters &parameters, const std::string &sys_name, const std::string &unknown_name), Gradient gptr(const Point &p, const Parameters &parameters, const std::string &sys_name, const std::string &unknown_name), Parameters &parameters, NumericVector< Number > &new_vector) const
unsigned int	number () const
void	update_global_solution (std::vector< Number > &global_soln) const
void	update_global_solution (std::vector< Number > &global_soln, const unsigned int dest_proc) const
const MeshBase &	get_mesh () const
MeshBase &	get_mesh ()
const DofMap &	get_dof_map () const
DofMap &	get_dof_map ()
const EquationSystems &	get_equation_systems () const
EquationSystems &	get_equation_systems ()
bool	active () const
void	activate ()
void	deactivate ()
vectors_iterator	vectors_begin ()
const_vectors_iterator	vectors_begin () const
vectors_iterator	vectors_end ()
const_vectors_iterator	vectors_end () const
NumericVector< Number > &	add_vector (const std::string &vec_name, const bool projections=true)
bool &	project_solution_on_reinit (void)
bool	have_vector (const std::string &vec_name) const
const NumericVector< Number > *	request_vector (const std::string &vec_name) const
NumericVector< Number > *	request_vector (const std::string &vec_name)
const NumericVector< Number > *	request_vector (const unsigned int vec_num) const
NumericVector< Number > *	request_vector (const unsigned int vec_num)
const NumericVector< Number > &	get_vector (const std::string &vec_name) const
NumericVector< Number > &	get_vector (const std::string &vec_name)
const NumericVector< Number > &	get_vector (const unsigned int vec_num) const
NumericVector< Number > &	get_vector (const unsigned int vec_num)
const std::string &	vector_name (const unsigned int vec_num)
NumericVector< Number > &	add_adjoint_solution (unsigned int i=0)
NumericVector< Number > &	get_adjoint_solution (unsigned int i=0)
const NumericVector< Number > &	get_adjoint_solution (unsigned int i=0) const
NumericVector< Number > &	add_sensitivity_solution (unsigned int i=0)
NumericVector< Number > &	get_sensitivity_solution (unsigned int i=0)
const NumericVector< Number > &	get_sensitivity_solution (unsigned int i=0) const
NumericVector< Number > &	add_weighted_sensitivity_adjoint_solution (unsigned int i=0)
NumericVector< Number > &	get_weighted_sensitivity_adjoint_solution (unsigned int i=0)
const NumericVector< Number > &	get_weighted_sensitivity_adjoint_solution (unsigned int i=0) const
NumericVector< Number > &	add_weighted_sensitivity_solution ()
NumericVector< Number > &	get_weighted_sensitivity_solution ()
const NumericVector< Number > &	get_weighted_sensitivity_solution () const
NumericVector< Number > &	add_adjoint_rhs (unsigned int i=0)
NumericVector< Number > &	get_adjoint_rhs (unsigned int i=0)
const NumericVector< Number > &	get_adjoint_rhs (unsigned int i=0) const
NumericVector< Number > &	add_sensitivity_rhs (unsigned int i=0)
NumericVector< Number > &	get_sensitivity_rhs (unsigned int i=0)
const NumericVector< Number > &	get_sensitivity_rhs (unsigned int i=0) const
unsigned int	n_vectors () const
unsigned int	n_vars () const
unsigned int	n_dofs () const
unsigned int	n_active_dofs () const
unsigned int	n_constrained_dofs () const
unsigned int	n_local_dofs () const
unsigned int	add_variable (const std::string &var, const FEType &type, const std::set< subdomain_id_type > *const active_subdomains=NULL)
unsigned int	add_variable (const std::string &var, const Order order=FIRST, const FEFamily=LAGRANGE, const std::set< subdomain_id_type > *const active_subdomains=NULL)
const Variable &	variable (unsigned int var) const
bool	has_variable (const std::string &var) const
const std::string &	variable_name (const unsigned int i) const
unsigned short int	variable_number (const std::string &var) const
const FEType &	variable_type (const unsigned int i) const
const FEType &	variable_type (const std::string &var) const
Real	calculate_norm (NumericVector< Number > &v, unsigned int var=0, FEMNormType norm_type=L2) const
Real	calculate_norm (NumericVector< Number > &v, const SystemNorm &norm) const
void	read_header (Xdr &io, const std::string &version, const bool read_header=true, const bool read_additional_data=true, const bool read_legacy_format=false)
void	read_legacy_data (Xdr &io, const bool read_additional_data=true)
void	read_serialized_data (Xdr &io, const bool read_additional_data=true)
void	read_parallel_data (Xdr &io, const bool read_additional_data)
void	write_header (Xdr &io, const std::string &version, const bool write_additional_data) const
void	write_serialized_data (Xdr &io, const bool write_additional_data=true) const
void	write_parallel_data (Xdr &io, const bool write_additional_data) const
std::string	get_info () const
void	attach_init_function (void fptr(EquationSystems &es, const std::string &name))
void	attach_assemble_function (void fptr(EquationSystems &es, const std::string &name))
void	attach_constraint_function (void fptr(EquationSystems &es, const std::string &name))
void	attach_QOI_function (void fptr(EquationSystems &es, const std::string &name, const QoISet &qoi_indices))
void	attach_QOI_derivative (void fptr(EquationSystems &es, const std::string &name, const QoISet &qoi_indices))
virtual void	user_initialization ()
virtual void	user_assembly ()
virtual void	user_constrain ()
virtual void	user_QOI (const QoISet &qoi_indices)
virtual void	user_QOI_derivative (const QoISet &qoi_indices)
virtual void	re_update ()
virtual void	restrict_vectors ()
virtual void	prolong_vectors ()
Number	current_solution (const unsigned int global_dof_number) const
void	local_dof_indices (const unsigned int var, std::set< unsigned int > &var_indices) const
void	zero_variable (NumericVector< Number > &v, unsigned int var_num) const
Static Public Member Functions
static std::string	get_info ()
static void	print_info ()
static unsigned int	n_objects ()
Public Attributes
AutoPtr< LinearSolver< Number > >	linear_solver
SparseMatrix< Number > *	matrix
NumericVector< Number > *	rhs
bool	assemble_before_solve
AutoPtr< NumericVector< Number > >	solution
AutoPtr< NumericVector< Number > >	current_local_solution
std::vector< Number >	qoi
Protected Types
typedef std::map< std::string, std::pair< unsigned int, unsigned int > >	Counts
Protected Member Functions
virtual void	init_data ()
virtual void	init_matrices ()
void	project_vector (NumericVector< Number > &) const
void	project_vector (const NumericVector< Number > &, NumericVector< Number > &) const
void	increment_constructor_count (const std::string &name)
void	increment_destructor_count (const std::string &name)
Protected Attributes
unsigned int	_n_linear_iterations
Real	_final_linear_residual
ShellMatrix< Number > *	_shell_matrix
Static Protected Attributes
static Counts	_counts
static Threads::atomic < unsigned int >	_n_objects
static Threads::spin_mutex	_mutex