oomph::FoldHandler Class Reference

#include <assembly_handler.h>

Inheritance diagram for oomph::FoldHandler:

Public Member Functions
	FoldHandler (Problem *const &problem_pt, double *const &parameter_pt)
	FoldHandler (Problem *const &problem_pt, double *const &parameter_pt, const DoubleVector &eigenvector)
	Constructor in which initial eigenvector can be passed. More...
	FoldHandler (Problem *const &problem_pt, double *const &parameter_pt, const DoubleVector &eigenvector, const DoubleVector &normalisation)
	~FoldHandler ()
	Destructor return the problem to its original (non-augmented) state. More...
unsigned	ndof (GeneralisedElement *const &elem_pt)
	Get the number of elemental degrees of freedom. More...
unsigned long	eqn_number (GeneralisedElement *const &elem_pt, const unsigned &ieqn_local)
	Get the global equation number of the local unknown. More...
void	get_residuals (GeneralisedElement *const &elem_pt, Vector< double > &residuals)
	Get the residuals. More...
void	get_jacobian (GeneralisedElement *const &elem_pt, Vector< double > &residuals, DenseMatrix< double > &jacobian)
void	get_dresiduals_dparameter (GeneralisedElement *const &elem_pt, double *const &parameter_pt, Vector< double > &dres_dparam)
void	get_djacobian_dparameter (GeneralisedElement *const &elem_pt, double *const &parameter_pt, Vector< double > &dres_dparam, DenseMatrix< double > &djac_dparam)
void	get_hessian_vector_products (GeneralisedElement *const &elem_pt, Vector< double > const &Y, DenseMatrix< double > const &C, DenseMatrix< double > &product)
int	bifurcation_type () const
	Indicate that we are tracking a fold bifurcation by returning 1. More...
double *	bifurcation_parameter_pt () const
void	get_eigenfunction (Vector< DoubleVector > &eigenfunction)
void	solve_augmented_block_system ()
	Set to solve the augmented block system. More...
void	solve_block_system ()
	Set to solve the block system. More...
void	solve_full_system ()
	Solve non-block system. More...
Public Member Functions inherited from oomph::AssemblyHandler
	AssemblyHandler ()
	Empty constructor. More...
virtual void	dof_vector (GeneralisedElement *const &elem_pt, const unsigned &t, Vector< double > &dof)
	Return vector of dofs at time level t in the element elem_pt. More...
virtual void	dof_pt_vector (GeneralisedElement *const &elem_pt, Vector< double * > &dof_pt)
	Return vector of pointers to dofs in the element elem_pt. More...
virtual double &	local_problem_dof (Problem *const &problem_pt, const unsigned &t, const unsigned &i)
virtual void	get_all_vectors_and_matrices (GeneralisedElement *const &elem_pt, Vector< Vector< double >> &vec, Vector< DenseMatrix< double >> &matrix)
virtual void	get_inner_products (GeneralisedElement *const &elem_pt, Vector< std::pair< unsigned, unsigned >> const &history_index, Vector< double > &inner_product)
virtual void	get_inner_product_vectors (GeneralisedElement *const &elem_pt, Vector< unsigned > const &history_index, Vector< Vector< double >> &inner_product_vector)
virtual	~AssemblyHandler ()
	Empty virtual destructor. More...

Private Types
enum	{ Full_augmented , Block_J , Block_augmented_J }

Private Attributes
unsigned	Solve_which_system
Problem *	Problem_pt
	Pointer to the problem. More...
unsigned	Ndof
Vector< double >	Phi
Vector< double >	Y
	Storage for the null vector. More...
Vector< int >	Count
double *	Parameter_pt
	Storage for the pointer to the parameter. More...

Friends
class	AugmentedBlockFoldLinearSolver

Detailed Description

A class that is used to assemble the augmented system that defines a fold (saddle-node) or limit point. The "standard" problem must be a function of a global paramter \(\lambda\), and a solution is \(R(u,\lambda) = 0 \) , where \( u \) are the unknowns in the problem. A limit point is formally specified by the augmented system of size \( 2N+1 \)

\[ R(u,\lambda) = 0, \]

\[ Jy = 0, \]

\[\phi\cdot y = 1. \]

In the above \( J \) is the usual Jacobian matrix, \( dR/du \), and \( \phi \) is a constant vector that is chosen to ensure that the null vector, \( y \), is not trivial.

Member Enumeration Documentation

anonymous enum

anonymous enum

private

A little private enum to determine whether we are solving the block system or not

Enumerator
Full_augmented
Block_J
Block_augmented_J

    {
      Full_augmented,
      Block_J,
      Block_augmented_J
    };

Constructor & Destructor Documentation

FoldHandler() [1/3]

oomph::FoldHandler::FoldHandler	(	Problem *const &	problem_pt,
		double *const &	parameter_pt
	)

Constructor: initialise the fold handler, by setting initial guesses for Y, Phi and calculating count. If the system changes, a new fold handler must be constructed

Constructor: Initialise the fold handler, by setting initial guesses for Y, Phi and calculating Count. If the system changes, a new handler must be constructed.

    : Solve_which_system(Full_augmented), Parameter_pt(parameter_pt)
  {
    // Set the problem pointer
    Problem_pt = problem_pt;
    // Set the number of degrees of freedom
    Ndof = problem_pt->ndof();
 
    // create the linear algebra distribution for this solver
    // currently only global (non-distributed) distributions are allowed
    LinearAlgebraDistribution* dist_pt =
      new LinearAlgebraDistribution(problem_pt->communicator_pt(), Ndof, false);
 
    // Resize the vectors of additional dofs and constants
    Phi.resize(Ndof);
    Y.resize(Ndof);
    Count.resize(Ndof, 0);
 
    // Loop over all the elements in the problem
    unsigned n_element = problem_pt->mesh_pt()->nelement();
    for (unsigned e = 0; e < n_element; e++)
    {
      GeneralisedElement* elem_pt = problem_pt->mesh_pt()->element_pt(e);
      // Loop over the local freedoms in the element
      unsigned n_var = elem_pt->ndof();
      for (unsigned n = 0; n < n_var; n++)
      {
        // Increase the associated global equation number counter
        ++Count[elem_pt->eqn_number(n)];
      }
    }
 
    // Calculate the value Phi by
    // solving the system JPhi = dF/dlambda
 
    // Locally cache the linear solver
    LinearSolver* const linear_solver_pt = problem_pt->linear_solver_pt();
 
    // Save the status before entry to this routine
    bool enable_resolve = linear_solver_pt->is_resolve_enabled();
 
    // We need to do a resolve
    linear_solver_pt->enable_resolve();
 
    // Storage for the solution
    DoubleVector x(dist_pt, 0.0);
 
    // Solve the standard problem, we only want to make sure that
    // we factorise the matrix, if it has not been factorised. We shall
    // ignore the return value of x.
    linear_solver_pt->solve(problem_pt, x);
 
    // Get the vector dresiduals/dparameter
    problem_pt->get_derivative_wrt_global_parameter(parameter_pt, x);
 
    // Copy rhs vector into local storage so it doesn't get overwritten
    // if the linear solver decides to initialise the solution vector, say,
    // which it's quite entitled to do!
    DoubleVector input_x(x);
 
    // Now resolve the system with the new RHS and overwrite the solution
    linear_solver_pt->resolve(input_x, x);
 
    // Restore the storage status of the linear solver
    if (enable_resolve)
    {
      linear_solver_pt->enable_resolve();
    }
    else
    {
      linear_solver_pt->disable_resolve();
    }
 
    // Add the global parameter as an unknown to the problem
    problem_pt->Dof_pt.push_back(parameter_pt);
 
 
    // Normalise the initial guesses for phi
    double length = 0.0;
    for (unsigned n = 0; n < Ndof; n++)
    {
      length += x[n] * x[n];
    }
    length = sqrt(length);
 
    // Now add the null space components to the problem unknowns
    // and initialise them and Phi to the same normalised values
    for (unsigned n = 0; n < Ndof; n++)
    {
      problem_pt->Dof_pt.push_back(&Y[n]);
      Y[n] = Phi[n] = -x[n] / length;
    }
 
    // delete the dist_pt
    problem_pt->Dof_distribution_pt->build(
      problem_pt->communicator_pt(), Ndof * 2 + 1, true);
    // Remove all previous sparse storage used during Jacobian assembly
    Problem_pt->Sparse_assemble_with_arrays_previous_allocation.resize(0);
 
    delete dist_pt;
  }

References oomph::LinearAlgebraDistribution::build(), oomph::Problem::communicator_pt(), Count, oomph::LinearSolver::disable_resolve(), oomph::Problem::Dof_distribution_pt, oomph::Problem::Dof_pt, e(), oomph::Mesh::element_pt(), oomph::LinearSolver::enable_resolve(), oomph::GeneralisedElement::eqn_number(), oomph::Problem::get_derivative_wrt_global_parameter(), oomph::LinearSolver::is_resolve_enabled(), oomph::Problem::linear_solver_pt(), oomph::Problem::mesh_pt(), n, Ndof, oomph::GeneralisedElement::ndof(), oomph::Problem::ndof(), oomph::Mesh::nelement(), Phi, Problem_pt, oomph::LinearSolver::resolve(), oomph::LinearSolver::solve(), oomph::Problem::Sparse_assemble_with_arrays_previous_allocation, sqrt(), plotDoE::x, and Y.

FoldHandler() [2/3]

oomph::FoldHandler::FoldHandler	(	Problem *const &	problem_pt,
		double *const &	parameter_pt,
		const DoubleVector &	eigenvector
	)

Constructor in which initial eigenvector can be passed.

Constructor in which the eigenvector is passed as an initial guess

    : Solve_which_system(Full_augmented), Parameter_pt(parameter_pt)
  {
    // Set the problem pointer
    Problem_pt = problem_pt;
    // Set the number of degrees of freedom
    Ndof = problem_pt->ndof();
 
    // create the linear algebra distribution for this solver
    // currently only global (non-distributed) distributions are allowed
    LinearAlgebraDistribution* dist_pt =
      new LinearAlgebraDistribution(problem_pt->communicator_pt(), Ndof, false);
 
    // Resize the vectors of additional dofs and constants
    Phi.resize(Ndof);
    Y.resize(Ndof);
    Count.resize(Ndof, 0);
 
    // Loop over all the elements in the problem
    unsigned n_element = problem_pt->mesh_pt()->nelement();
    for (unsigned e = 0; e < n_element; e++)
    {
      GeneralisedElement* elem_pt = problem_pt->mesh_pt()->element_pt(e);
      // Loop over the local freedoms in the element
      unsigned n_var = elem_pt->ndof();
      for (unsigned n = 0; n < n_var; n++)
      {
        // Increase the associated global equation number counter
        ++Count[elem_pt->eqn_number(n)];
      }
    }
 
    // Add the global parameter as an unknown to the problem
    problem_pt->Dof_pt.push_back(parameter_pt);
 
 
    // Normalise the initial guesses for the eigenvecor
    double length = 0.0;
    for (unsigned n = 0; n < Ndof; n++)
    {
      length += eigenvector[n] * eigenvector[n];
    }
    length = sqrt(length);
 
    // Now add the null space components to the problem unknowns
    // and initialise them and Phi to the same normalised values
    for (unsigned n = 0; n < Ndof; n++)
    {
      problem_pt->Dof_pt.push_back(&Y[n]);
      Y[n] = Phi[n] = eigenvector[n] / length;
    }
 
    // delete the dist_pt
    problem_pt->Dof_distribution_pt->build(
      problem_pt->communicator_pt(), Ndof * 2 + 1, true);
    // Remove all previous sparse storage used during Jacobian assembly
    Problem_pt->Sparse_assemble_with_arrays_previous_allocation.resize(0);
 
    delete dist_pt;
  }

References oomph::LinearAlgebraDistribution::build(), oomph::Problem::communicator_pt(), Count, oomph::Problem::Dof_distribution_pt, oomph::Problem::Dof_pt, e(), oomph::Mesh::element_pt(), oomph::GeneralisedElement::eqn_number(), oomph::Problem::mesh_pt(), n, Ndof, oomph::GeneralisedElement::ndof(), oomph::Problem::ndof(), oomph::Mesh::nelement(), Phi, Problem_pt, oomph::Problem::Sparse_assemble_with_arrays_previous_allocation, sqrt(), and Y.

FoldHandler() [3/3]

oomph::FoldHandler::FoldHandler	(	Problem *const &	problem_pt,
		double *const &	parameter_pt,
		const DoubleVector &	eigenvector,
		const DoubleVector &	normalisation
	)

Constructor in which initial eigenvector and normalisation can be passed

Constructor in which the eigenvector and normalisation vector are passed as an initial guess

    : Solve_which_system(Full_augmented), Parameter_pt(parameter_pt)
  {
    // Set the problem pointer
    Problem_pt = problem_pt;
    // Set the number of degrees of freedom
    Ndof = problem_pt->ndof();
 
    // create the linear algebra distribution for this solver
    // currently only global (non-distributed) distributions are allowed
    LinearAlgebraDistribution* dist_pt =
      new LinearAlgebraDistribution(problem_pt->communicator_pt(), Ndof, false);
 
    // Resize the vectors of additional dofs and constants
    Phi.resize(Ndof);
    Y.resize(Ndof);
    Count.resize(Ndof, 0);
 
    // Loop over all the elements in the problem
    unsigned n_element = problem_pt->mesh_pt()->nelement();
    for (unsigned e = 0; e < n_element; e++)
    {
      GeneralisedElement* elem_pt = problem_pt->mesh_pt()->element_pt(e);
      // Loop over the local freedoms in the element
      unsigned n_var = elem_pt->ndof();
      for (unsigned n = 0; n < n_var; n++)
      {
        // Increase the associated global equation number counter
        ++Count[elem_pt->eqn_number(n)];
      }
    }
 
    // Add the global parameter as an unknown to the problem
    problem_pt->Dof_pt.push_back(parameter_pt);
 
 
    // Normalise the initial guesses for the eigenvecor
    double length = 0.0;
    for (unsigned n = 0; n < Ndof; n++)
    {
      length += eigenvector[n] * normalisation[n];
    }
    length = sqrt(length);
 
    // Now add the null space components to the problem unknowns
    // and initialise them and Phi to the same normalised values
    for (unsigned n = 0; n < Ndof; n++)
    {
      problem_pt->Dof_pt.push_back(&Y[n]);
      Y[n] = eigenvector[n] / length;
      Phi[n] = normalisation[n];
    }
 
    // delete the dist_pt
    problem_pt->Dof_distribution_pt->build(
      problem_pt->communicator_pt(), Ndof * 2 + 1, true);
    // Remove all previous sparse storage used during Jacobian assembly
    Problem_pt->Sparse_assemble_with_arrays_previous_allocation.resize(0);
 
    delete dist_pt;
  }

References oomph::LinearAlgebraDistribution::build(), oomph::Problem::communicator_pt(), Count, oomph::Problem::Dof_distribution_pt, oomph::Problem::Dof_pt, e(), oomph::Mesh::element_pt(), oomph::GeneralisedElement::eqn_number(), oomph::Problem::mesh_pt(), n, Ndof, oomph::GeneralisedElement::ndof(), oomph::Problem::ndof(), oomph::Mesh::nelement(), Phi, Problem_pt, oomph::Problem::Sparse_assemble_with_arrays_previous_allocation, sqrt(), and Y.

~FoldHandler()

oomph::FoldHandler::~FoldHandler

(

)

Destructor return the problem to its original (non-augmented) state.

Destructor, return the problem to its original state before the augmented system was added

  {
    // If we are using the block solver reset the problem's linear solver
    // to the original one
    AugmentedBlockFoldLinearSolver* block_fold_solver_pt =
      dynamic_cast<AugmentedBlockFoldLinearSolver*>(
        Problem_pt->linear_solver_pt());
 
    if (block_fold_solver_pt)
    {
      // Reset the problem's linear solver
      Problem_pt->linear_solver_pt() = block_fold_solver_pt->linear_solver_pt();
      // Delete the block solver
      delete block_fold_solver_pt;
    }
 
    // Resize the number of dofs
    Problem_pt->Dof_pt.resize(Ndof);
    Problem_pt->Dof_distribution_pt->build(
      Problem_pt->communicator_pt(), Ndof, false);
    // Remove all previous sparse storage used during Jacobian assembly
    Problem_pt->Sparse_assemble_with_arrays_previous_allocation.resize(0);
  }

References oomph::LinearAlgebraDistribution::build(), oomph::Problem::communicator_pt(), oomph::Problem::Dof_distribution_pt, oomph::Problem::Dof_pt, oomph::Problem::linear_solver_pt(), oomph::AugmentedBlockFoldLinearSolver::linear_solver_pt(), Ndof, Problem_pt, and oomph::Problem::Sparse_assemble_with_arrays_previous_allocation.

Member Function Documentation

bifurcation_parameter_pt()

double* oomph::FoldHandler::bifurcation_parameter_pt ( ) const

inlinevirtual

Return a pointer to the bifurcation parameter in bifurcation tracking problems

Reimplemented from oomph::AssemblyHandler.

    {
      return Parameter_pt;
    }

References Parameter_pt.

bifurcation_type()

int oomph::FoldHandler::bifurcation_type ( ) const

inlinevirtual

Indicate that we are tracking a fold bifurcation by returning 1.

Reimplemented from oomph::AssemblyHandler.

    {
      return 1;
    }

eqn_number()

unsigned long oomph::FoldHandler::eqn_number ( GeneralisedElement *const &  elem_pt, const unsigned &  ieqn_local  )

virtual

Get the global equation number of the local unknown.

Return the global equation number associated with local equation number ieqn_local. We choose to number the unknowns according to the augmented system.

Reimplemented from oomph::AssemblyHandler.

  {
    // Find the number of non-augmented dofs in the element
    unsigned raw_ndof = elem_pt->ndof();
    // Storage for the global eqn number
    unsigned long global_eqn = 0;
    // If we are a "standard" unknown, just return the standard equation number
    if (ieqn_local < raw_ndof)
    {
      global_eqn = elem_pt->eqn_number(ieqn_local);
    }
    // Otherwise if we are at an unknown corresponding to the bifurcation
    // parameter return the global equation number of the parameter
    else if (ieqn_local == raw_ndof)
    {
      global_eqn = Ndof;
    }
    // Otherwise we are in the unknown corresponding to a null vector
    // return the global unknown Ndof + 1 + global unknown of "standard"
    // unknown.
    else
    {
      global_eqn = Ndof + 1 + elem_pt->eqn_number(ieqn_local - 1 - raw_ndof);
    }
 
    // Return the global equation number
    return global_eqn;
  }

References oomph::GeneralisedElement::eqn_number(), Ndof, and oomph::GeneralisedElement::ndof().

Referenced by get_jacobian(), oomph::AugmentedBlockFoldLinearSolver::resolve(), and oomph::AugmentedBlockFoldLinearSolver::solve().

get_djacobian_dparameter()

void oomph::FoldHandler::get_djacobian_dparameter ( GeneralisedElement *const &  elem_pt, double *const &  parameter_pt, Vector< double > &  dres_dparam, DenseMatrix< double > &  djac_dparam  )

virtual

Overload the derivative of the residuals and jacobian with respect to a parameter so that it breaks

Overload the derivative of the residuals and jacobian with respect to a parameter so that it breaks because it should not be required

Reimplemented from oomph::AssemblyHandler.

  {
    std::ostringstream error_stream;
    error_stream
      << "This function has not been implemented because it is not required\n";
    error_stream << "in standard problems.\n";
    error_stream
      << "If you find that you need it, you will have to implement it!\n\n";
 
    throw OomphLibError(
      error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
  }

References OOMPH_CURRENT_FUNCTION, and OOMPH_EXCEPTION_LOCATION.

get_dresiduals_dparameter()

void oomph::FoldHandler::get_dresiduals_dparameter ( GeneralisedElement *const &  elem_pt, double *const &  parameter_pt, Vector< double > &  dres_dparam  )

virtual

Overload the derivatives of the residuals with respect to a parameter to apply to the augmented system

Formulate the derivatives of the augmented system with respect to a parameter

Reimplemented from oomph::AssemblyHandler.

  {
    // Need to get raw residuals and jacobian
    unsigned raw_ndof = elem_pt->ndof();
 
    // Find out which system we are solving
    switch (Solve_which_system)
    {
        // If we are solving the standard system
      case Block_J:
      {
        // Get the basic residual derivatives
        elem_pt->get_dresiduals_dparameter(parameter_pt, dres_dparam);
      }
      break;
 
      // If we are solving the augmented-by-one system
      case Block_augmented_J:
      {
        // Get the basic residual derivatives
        elem_pt->get_dresiduals_dparameter(parameter_pt, dres_dparam);
 
        // Zero the final derivative
        dres_dparam[raw_ndof] = 0.0;
      }
      break;
 
      // If we are solving the full augmented system
      case Full_augmented:
      {
        DenseMatrix<double> djac_dparam(raw_ndof);
        // Get the basic residuals and jacobian derivatives initially
        elem_pt->get_djacobian_dparameter(
          parameter_pt, dres_dparam, djac_dparam);
 
        // The normalisation equation does not depend on the parameter
        dres_dparam[raw_ndof] = 0.0;
 
        // Now assemble the equations dJy/dparameter = 0
        for (unsigned i = 0; i < raw_ndof; i++)
        {
          dres_dparam[raw_ndof + 1 + i] = 0.0;
          for (unsigned j = 0; j < raw_ndof; j++)
          {
            dres_dparam[raw_ndof + 1 + i] +=
              djac_dparam(i, j) * Y[elem_pt->eqn_number(j)];
          }
        }
      }
      break;
 
      default:
        std::ostringstream error_stream;
        error_stream
          << "The Solve_which_system flag can only take values 0, 1, 2"
          << " not " << Solve_which_system << "\n";
        throw OomphLibError(
          error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
    }
  }

References Block_augmented_J, Block_J, oomph::GeneralisedElement::eqn_number(), Full_augmented, oomph::GeneralisedElement::get_djacobian_dparameter(), oomph::GeneralisedElement::get_dresiduals_dparameter(), i, j, oomph::GeneralisedElement::ndof(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, Solve_which_system, and Y.

get_eigenfunction()

void oomph::FoldHandler::get_eigenfunction ( Vector< DoubleVector > &  eigenfunction )

virtual

Return the eigenfunction(s) associated with the bifurcation that has been detected in bifurcation tracking problems

Reimplemented from oomph::AssemblyHandler.

  {
    // There is only one (real) null vector
    eigenfunction.resize(1);
    LinearAlgebraDistribution dist(Problem_pt->communicator_pt(), Ndof, false);
    // Rebuild the vector
    eigenfunction[0].build(&dist, 0.0);
    // Set the value to be the null vector
    for (unsigned n = 0; n < Ndof; n++)
    {
      eigenfunction[0][n] = Y[n];
    }
  }

References oomph::Problem::communicator_pt(), n, Ndof, Problem_pt, and Y.

get_hessian_vector_products()

void oomph::FoldHandler::get_hessian_vector_products ( GeneralisedElement *const &  elem_pt, Vector< double > const &  Y, DenseMatrix< double > const &  C, DenseMatrix< double > &  product  )

virtual

Overload the hessian vector product function so that it breaks

Overload the hessian vector product function so that it breaks because it should not be required

Reimplemented from oomph::AssemblyHandler.

  {
    std::ostringstream error_stream;
    error_stream
      << "This function has not been implemented because it is not required\n";
    error_stream << "in standard problems.\n";
    error_stream
      << "If you find that you need it, you will have to implement it!\n\n";
 
    throw OomphLibError(
      error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
  }

References OOMPH_CURRENT_FUNCTION, and OOMPH_EXCEPTION_LOCATION.

get_jacobian()

void oomph::FoldHandler::get_jacobian ( GeneralisedElement *const &  elem_pt, Vector< double > &  residuals, DenseMatrix< double > &  jacobian  )

virtual

Calculate the elemental Jacobian matrix "d equation / d variable".

Calculate the elemental Jacobian matrix "d equation / d variable" for the augmented system

ALICE

ALICE

Reimplemented from oomph::AssemblyHandler.

  {
    // Find the number of augmented dofs
    unsigned augmented_ndof = ndof(elem_pt);
    // Find the non-augmented dofs
    unsigned raw_ndof = elem_pt->ndof();
 
    // Which system are we solving
    switch (Solve_which_system)
    {
        // If we are solving the original system
      case Block_J:
      {
        // Just get the raw jacobian and residuals
        elem_pt->get_jacobian(residuals, jacobian);
      }
      break;
 
      // If we are solving the augmented-by-one system
      case Block_augmented_J:
      {
        // Get the full residuals, we need them
        get_residuals(elem_pt, residuals);
 
        // Need to get the raw jacobian (and raw residuals)
        Vector<double> newres(augmented_ndof);
        elem_pt->get_jacobian(newres, jacobian);
 
        // Now do finite differencing stuff
        const double FD_step = 1.0e-8;
        // Fill in the first lot of finite differences
        {
          // increase the global parameter
          double* unknown_pt = Problem_pt->Dof_pt[Ndof];
          double init = *unknown_pt;
          *unknown_pt += FD_step;
 
          // Now do any possible updates
          Problem_pt->actions_after_change_in_bifurcation_parameter();
 
          // Get the new (modified) residuals
          get_residuals(elem_pt, newres);
 
          // The final column  is given by the difference
          // between the residuals
          for (unsigned n = 0; n < raw_ndof; n++)
          {
            jacobian(n, augmented_ndof - 1) =
              (newres[n] - residuals[n]) / FD_step;
          }
          // Reset the global parameter
          *unknown_pt = init;
 
          // Now do any possible updates
          Problem_pt->actions_after_change_in_bifurcation_parameter();
        }
 
        // Fill in the bottom row
        for (unsigned n = 0; n < raw_ndof; n++)
        {
          unsigned local_eqn = elem_pt->eqn_number(n);
          jacobian(augmented_ndof - 1, n) = Phi[local_eqn] / Count[local_eqn];
        }
      }
      break;
 
        // Otherwise solving the full system
      case Full_augmented:
      {
        // Get the residuals for the augmented system
        get_residuals(elem_pt, residuals);
 
        // Need to get the raw residuals and jacobian
        Vector<double> newres(raw_ndof);
        DenseMatrix<double> newjac(raw_ndof);
        elem_pt->get_jacobian(newres, jacobian);
 
        // Fill in the jacobian on the diagonal sub-block of
        // the null-space equations
        for (unsigned n = 0; n < raw_ndof; n++)
        {
          for (unsigned m = 0; m < raw_ndof; m++)
          {
            jacobian(raw_ndof + 1 + n, raw_ndof + 1 + m) = jacobian(n, m);
          }
        }
 
        // Now finite difference wrt the global unknown
        const double FD_step = 1.0e-8;
        // Fill in the first lot of finite differences
        {
          // increase the global parameter
          double* unknown_pt = Problem_pt->Dof_pt[Ndof];
          double init = *unknown_pt;
          *unknown_pt += FD_step;
          // Need to update the function
          Problem_pt->actions_after_change_in_bifurcation_parameter();
 
          // Get the new raw residuals and jacobian
          elem_pt->get_jacobian(newres, newjac);
 
          // The end of the first row is given by the difference
          // between the residuals
          for (unsigned n = 0; n < raw_ndof; n++)
          {
            jacobian(n, raw_ndof) = (newres[n] - residuals[n]) / FD_step;
            // The end of the second row is given by the difference multiplied
            // by the product
            for (unsigned l = 0; l < raw_ndof; l++)
            {
              jacobian(raw_ndof + 1 + n, raw_ndof) +=
                (newjac(n, l) - jacobian(n, l)) * Y[elem_pt->eqn_number(l)] /
                FD_step;
            }
          }
          // Reset the global parameter
          *unknown_pt = init;
 
          // Need to update the function
          Problem_pt->actions_after_change_in_bifurcation_parameter();
        }
 
        // Now fill in the first column of the second rows
        {
          for (unsigned n = 0; n < raw_ndof; n++)
          {
            unsigned long global_eqn = eqn_number(elem_pt, n);
            // Increase the first lot
            double* unknown_pt = Problem_pt->Dof_pt[global_eqn];
            double init = *unknown_pt;
            *unknown_pt += FD_step;
            Problem_pt->actions_before_newton_convergence_check(); 
 
            // Get the new jacobian
            elem_pt->get_jacobian(newres, newjac);
 
            // Work out the differences
            for (unsigned k = 0; k < raw_ndof; k++)
            {
              // jacobian(raw_ndof+k,n) = 0.0;
              for (unsigned l = 0; l < raw_ndof; l++)
              {
                jacobian(raw_ndof + 1 + k, n) +=
                  (newjac(k, l) - jacobian(k, l)) * Y[elem_pt->eqn_number(l)] /
                  FD_step;
              }
            }
            *unknown_pt = init;
            Problem_pt->actions_before_newton_convergence_check(); 
          }
        }
 
        // Fill in the row corresponding to the parameter
        for (unsigned n = 0; n < raw_ndof; n++)
        {
          unsigned global_eqn = elem_pt->eqn_number(n);
          jacobian(raw_ndof, raw_ndof + 1 + n) =
            Phi[global_eqn] / Count[global_eqn];
        }
 
        //   //Loop over the dofs
        //     for(unsigned n=0;n<augmented_ndof;n++)
        //      {
        //       unsigned long global_eqn = eqn_number(elem_pt,n);
        //       double* unknown_pt = Problem_pt->Dof_pt[global_eqn];
        //       double init = *unknown_pt;
        //       *unknown_pt += FD_step;
 
        //       //Get the new residuals
        //       get_residuals(elem_pt,newres);
 
        //       for(unsigned m=0;m<augmented_ndof;m++)
        //        {
        //         jacobian(m,n) = (newres[m] - residuals[m])/FD_step;
        //        }
        //       //Reset the unknown
        //       *unknown_pt = init;
        //      }
      }
      break;
 
      default:
        std::ostringstream error_stream;
        error_stream
          << "The Solve_which_system flag can only take values 0, 1, 2"
          << " not " << Solve_which_system << "\n";
        throw OomphLibError(
          error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
    }
  }

References oomph::Problem::actions_after_change_in_bifurcation_parameter(), oomph::Problem::actions_before_newton_convergence_check(), Block_augmented_J, Block_J, Count, oomph::Problem::Dof_pt, oomph::GeneralisedElement::eqn_number(), eqn_number(), oomph::BlackBoxFDNewtonSolver::FD_step, Full_augmented, oomph::GeneralisedElement::get_jacobian(), get_residuals(), k, m, n, Ndof, oomph::GeneralisedElement::ndof(), ndof(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, Phi, Problem_pt, Solve_which_system, and Y.

Referenced by oomph::AugmentedBlockFoldLinearSolver::resolve(), and oomph::AugmentedBlockFoldLinearSolver::solve().

get_residuals()

void oomph::FoldHandler::get_residuals ( GeneralisedElement *const &  elem_pt, Vector< double > &  residuals  )

virtual

Get the residuals.

Formulate the augmented system.

Reimplemented from oomph::AssemblyHandler.

  {
    // Need to get raw residuals and jacobian
    unsigned raw_ndof = elem_pt->ndof();
 
    // Find out which system we are solving
    switch (Solve_which_system)
    {
        // If we are solving the standard system
      case Block_J:
      {
        // Get the basic residuals
        elem_pt->get_residuals(residuals);
      }
      break;
 
      // If we are solving the augmented-by-one system
      case Block_augmented_J:
      {
        // Get the basic residuals
        elem_pt->get_residuals(residuals);
 
        // Zero the final residual
        residuals[raw_ndof] = 0.0;
      }
      break;
 
      // If we are solving the full augmented system
      case Full_augmented:
      {
        DenseMatrix<double> jacobian(raw_ndof);
        // Get the basic residuals and jacobian initially
        elem_pt->get_jacobian(residuals, jacobian);
 
        // The normalisation equation must be initialised to -1.0/number of
        // elements
        residuals[raw_ndof] = -1.0 / Problem_pt->mesh_pt()->nelement();
 
        // Now assemble the equations Jy = 0
        for (unsigned i = 0; i < raw_ndof; i++)
        {
          residuals[raw_ndof + 1 + i] = 0.0;
          for (unsigned j = 0; j < raw_ndof; j++)
          {
            residuals[raw_ndof + 1 + i] +=
              jacobian(i, j) * Y[elem_pt->eqn_number(j)];
          }
          // Add the contribution to phi.y=1
          // Need to divide by the number of elements that contribute to this
          // unknown so that we do get phi.y exactly.
          unsigned global_eqn = elem_pt->eqn_number(i);
          residuals[raw_ndof] +=
            (Phi[global_eqn] * Y[global_eqn]) / Count[global_eqn];
        }
      }
      break;
 
      default:
        std::ostringstream error_stream;
        error_stream
          << "The Solve_which_system flag can only take values 0, 1, 2"
          << " not " << Solve_which_system << "\n";
        throw OomphLibError(
          error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
    }
  }

References Block_augmented_J, Block_J, Count, oomph::GeneralisedElement::eqn_number(), Full_augmented, oomph::GeneralisedElement::get_jacobian(), oomph::GeneralisedElement::get_residuals(), i, j, oomph::Problem::mesh_pt(), oomph::GeneralisedElement::ndof(), oomph::Mesh::nelement(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, Phi, Problem_pt, Solve_which_system, and Y.

Referenced by get_jacobian().

ndof()

unsigned oomph::FoldHandler::ndof ( GeneralisedElement *const &  elem_pt )

virtual

Get the number of elemental degrees of freedom.

The number of unknowns is 2n+1.

Reimplemented from oomph::AssemblyHandler.

  {
    unsigned raw_ndof = elem_pt->ndof();
    // Return different values depending on the type of block decomposition
    switch (Solve_which_system)
    {
      case Full_augmented:
        return (2 * raw_ndof + 1);
        break;
 
      case Block_augmented_J:
        return (raw_ndof + 1);
        break;
 
      case Block_J:
        return raw_ndof;
        break;
 
      default:
        std::ostringstream error_stream;
        error_stream
          << "The Solve_which_system flag can only take values 0, 1, 2"
          << " not " << Solve_which_system << "\n";
        throw OomphLibError(
          error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
    }
  }

References Block_augmented_J, Block_J, Full_augmented, oomph::GeneralisedElement::ndof(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, and Solve_which_system.

Referenced by get_jacobian(), oomph::AugmentedBlockFoldLinearSolver::resolve(), and oomph::AugmentedBlockFoldLinearSolver::solve().

solve_augmented_block_system()

void oomph::FoldHandler::solve_augmented_block_system

(

)

Set to solve the augmented block system.

Set to solve the augmented-by-one block system.

  {
    // Only bother to do anything if we haven't already set the flag
    if (Solve_which_system != Block_augmented_J)
    {
      // If we were solving the system with the original jacobian add the
      // parameter
      if (Solve_which_system == Block_J)
      {
        Problem_pt->Dof_pt.push_back(Parameter_pt);
      }
 
      // Restrict the problem to the standard variables and
      // the bifurcation parameter only
      Problem_pt->Dof_pt.resize(Ndof + 1);
      Problem_pt->Dof_distribution_pt->build(
        Problem_pt->communicator_pt(), Ndof + 1, false);
      // Remove all previous sparse storage used during Jacobian assembly
      Problem_pt->Sparse_assemble_with_arrays_previous_allocation.resize(0);
      // Set the solve flag
      Solve_which_system = Block_augmented_J;
    }
  }

References Block_augmented_J, Block_J, oomph::LinearAlgebraDistribution::build(), oomph::Problem::communicator_pt(), oomph::Problem::Dof_distribution_pt, oomph::Problem::Dof_pt, Ndof, Parameter_pt, Problem_pt, Solve_which_system, and oomph::Problem::Sparse_assemble_with_arrays_previous_allocation.

Referenced by oomph::AugmentedBlockFoldLinearSolver::resolve(), and oomph::AugmentedBlockFoldLinearSolver::solve().

solve_block_system()

void oomph::FoldHandler::solve_block_system

(

)

Set to solve the block system.

Set to solve the block system. The boolean flag specifies whether the block decomposition should use exactly the same jacobian

  {
    // Only bother to do anything if we haven't already set the flag
    if (Solve_which_system != Block_J)
    {
      // Restrict the problem to the standard variables
      Problem_pt->Dof_pt.resize(Ndof);
      Problem_pt->Dof_distribution_pt->build(
        Problem_pt->communicator_pt(), Ndof, false);
      // Remove all previous sparse storage used during Jacobian assembly
      Problem_pt->Sparse_assemble_with_arrays_previous_allocation.resize(0);
 
      // Set the solve flag
      Solve_which_system = Block_J;
    }
  }

References Block_J, oomph::LinearAlgebraDistribution::build(), oomph::Problem::communicator_pt(), oomph::Problem::Dof_distribution_pt, oomph::Problem::Dof_pt, Ndof, Problem_pt, Solve_which_system, and oomph::Problem::Sparse_assemble_with_arrays_previous_allocation.

solve_full_system()

void oomph::FoldHandler::solve_full_system

(

)

Solve non-block system.

Set to Solve non-block system.

  {
    // Only do something if we are not solving the full system
    if (Solve_which_system != Full_augmented)
    {
      // If we were solving the problem without any augmentation,
      // add the parameter again
      if (Solve_which_system == Block_J)
      {
        Problem_pt->Dof_pt.push_back(Parameter_pt);
      }
 
      // Always add the additional unknowns back into the problem
      for (unsigned n = 0; n < Ndof; n++)
      {
        Problem_pt->Dof_pt.push_back(&Y[n]);
      }
 
      // update the Dof distribution pt
      Problem_pt->Dof_distribution_pt->build(
        Problem_pt->communicator_pt(), Ndof * 2 + 1, false);
      // Remove all previous sparse storage used during Jacobian assembly
      Problem_pt->Sparse_assemble_with_arrays_previous_allocation.resize(0);
 
      Solve_which_system = Full_augmented;
    }
  }

References Block_J, oomph::LinearAlgebraDistribution::build(), oomph::Problem::communicator_pt(), oomph::Problem::Dof_distribution_pt, oomph::Problem::Dof_pt, Full_augmented, n, Ndof, Parameter_pt, Problem_pt, Solve_which_system, oomph::Problem::Sparse_assemble_with_arrays_previous_allocation, and Y.

Referenced by oomph::AugmentedBlockFoldLinearSolver::resolve(), and oomph::AugmentedBlockFoldLinearSolver::solve().

Friends And Related Function Documentation

AugmentedBlockFoldLinearSolver

friend class AugmentedBlockFoldLinearSolver

friend

Member Data Documentation

Count

Vector<int> oomph::FoldHandler::Count

private

A vector that is used to determine how many elements contribute to a particular equation. It is used to ensure that the global system is correctly formulated.

Referenced by FoldHandler(), get_jacobian(), and get_residuals().

Ndof

unsigned oomph::FoldHandler::Ndof

private

Store the number of degrees of freedom in the non-augmented problem

Referenced by eqn_number(), FoldHandler(), get_eigenfunction(), get_jacobian(), solve_augmented_block_system(), solve_block_system(), solve_full_system(), and ~FoldHandler().

Parameter_pt

double* oomph::FoldHandler::Parameter_pt

private

Storage for the pointer to the parameter.

Referenced by bifurcation_parameter_pt(), solve_augmented_block_system(), and solve_full_system().

Phi

Vector<double> oomph::FoldHandler::Phi

private

A constant vector used to ensure that the null vector is not trivial

Referenced by FoldHandler(), get_jacobian(), and get_residuals().

Problem_pt

Problem* oomph::FoldHandler::Problem_pt

private

Pointer to the problem.

Referenced by FoldHandler(), get_eigenfunction(), get_jacobian(), get_residuals(), solve_augmented_block_system(), solve_block_system(), solve_full_system(), and ~FoldHandler().

Solve_which_system

unsigned oomph::FoldHandler::Solve_which_system

private

Integer flag to indicate which system should be assembled. There are three possibilities. The full augmented system (0), the non-augmented jacobian system (1), a system in which the jacobian is augmented by 1 row and column to ensure that it is non-singular (2). See the enum above

Referenced by get_dresiduals_dparameter(), get_jacobian(), get_residuals(), ndof(), solve_augmented_block_system(), solve_block_system(), and solve_full_system().

Y

Vector<double> oomph::FoldHandler::Y

private

Storage for the null vector.

Referenced by FoldHandler(), get_dresiduals_dparameter(), get_eigenfunction(), get_jacobian(), get_residuals(), oomph::AugmentedBlockFoldLinearSolver::resolve(), oomph::AugmentedBlockFoldLinearSolver::solve(), and solve_full_system().

---

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

• assembly_handler.h
• assembly_handler.cc