assembly_handler.h

Go to the documentation of this file.

// LIC// ====================================================================
// LIC// This file forms part of oomph-lib, the object-oriented,
// LIC// multi-physics finite-element library, available
// LIC// at http://www.oomph-lib.org.
// LIC//
// LIC// Copyright (C) 2006-2022 Matthias Heil and Andrew Hazel
// LIC//
// LIC// This library is free software; you can redistribute it and/or
// LIC// modify it under the terms of the GNU Lesser General Public
// LIC// License as published by the Free Software Foundation; either
// LIC// version 2.1 of the License, or (at your option) any later version.
// LIC//
// LIC// This library is distributed in the hope that it will be useful,
// LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
// LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
// LIC// Lesser General Public License for more details.
// LIC//
// LIC// You should have received a copy of the GNU Lesser General Public
// LIC// License along with this library; if not, write to the Free Software
// LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
// LIC// 02110-1301  USA.
// LIC//
// LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
// LIC//
// LIC//====================================================================
// A class that is used to assemble the linear systems solved
// by oomph-lib.
 
// Include guards to prevent multiple inclusion of this header
#ifndef OOMPH_ASSEMBLY_HANDLER_CLASS_HEADER
#define OOMPH_ASSEMBLY_HANDLER_CLASS_HEADER
 
// Config header generated by autoconfig
#ifdef HAVE_CONFIG_H
#include <oomph-lib-config.h>
#endif
 
// OOMPH-LIB headers
#include "matrices.h"
#include "linear_solver.h"
#include "double_vector_with_halo.h"
 
namespace oomph
{
  // Forward class definition of the element
  class GeneralisedElement;
 
  // Forward class definition of the problem
  class Problem;
 
  //=============================================================
  //===============================================================
  class AssemblyHandler
  {
  public:
    AssemblyHandler() {}
 
    virtual unsigned ndof(GeneralisedElement* const& elem_pt);
 
    virtual void dof_vector(GeneralisedElement* const& elem_pt,
                            const unsigned& t,
                            Vector<double>& dof);
 
    virtual void dof_pt_vector(GeneralisedElement* const& elem_pt,
                               Vector<double*>& dof_pt);
 
    virtual double& local_problem_dof(Problem* const& problem_pt,
                                      const unsigned& t,
                                      const unsigned& i);
 
    virtual unsigned long eqn_number(GeneralisedElement* const& elem_pt,
                                     const unsigned& ieqn_local);
 
    virtual void get_residuals(GeneralisedElement* const& elem_pt,
                               Vector<double>& residuals);
 
    virtual void get_jacobian(GeneralisedElement* const& elem_pt,
                              Vector<double>& residuals,
                              DenseMatrix<double>& jacobian);
 
    virtual void get_all_vectors_and_matrices(
      GeneralisedElement* const& elem_pt,
      Vector<Vector<double>>& vec,
      Vector<DenseMatrix<double>>& matrix);
 
    virtual void get_dresiduals_dparameter(GeneralisedElement* const& elem_pt,
                                           double* const& parameter_pt,
                                           Vector<double>& dres_dparam);
 
    virtual void get_djacobian_dparameter(GeneralisedElement* const& elem_pt,
                                          double* const& parameter_pt,
                                          Vector<double>& dres_dparam,
                                          DenseMatrix<double>& djac_dparam);
 
    virtual void get_hessian_vector_products(GeneralisedElement* const& elem_pt,
                                             Vector<double> const& Y,
                                             DenseMatrix<double> const& C,
                                             DenseMatrix<double>& product);
 
 
    virtual int bifurcation_type() const
    {
      return 0;
    }
 
    virtual double* bifurcation_parameter_pt() const;
 
    virtual void get_eigenfunction(Vector<DoubleVector>& eigenfunction);
 
    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);
 
#ifdef OOMPH_HAS_MPI
 
    virtual void synchronise() {}
#endif
 
    virtual ~AssemblyHandler() {}
  };
 
 
  //=============================================================
  //===============================================================
  class ExplicitTimeStepHandler : public AssemblyHandler
  {
  public:
    ExplicitTimeStepHandler() {}
 
    unsigned ndof(GeneralisedElement* const& elem_pt);
 
    unsigned long eqn_number(GeneralisedElement* const& elem_pt,
                             const unsigned& ieqn_local);
 
    void get_residuals(GeneralisedElement* const& elem_pt,
                       Vector<double>& residuals);
 
    void get_jacobian(GeneralisedElement* const& elem_pt,
                      Vector<double>& residuals,
                      DenseMatrix<double>& jacobian);
 
    void get_all_vectors_and_matrices(GeneralisedElement* const& elem_pt,
                                      Vector<Vector<double>>& vec,
                                      Vector<DenseMatrix<double>>& matrix);
 
    ~ExplicitTimeStepHandler() {}
  };
 
 
  //=============================================================
  //===============================================================
  class EigenProblemHandler : public AssemblyHandler
  {
    double Sigma_real;
 
  public:
    EigenProblemHandler(const double& sigma_real) : Sigma_real(sigma_real) {}
 
    unsigned ndof(GeneralisedElement* const& elem_pt);
 
    unsigned long eqn_number(GeneralisedElement* const& elem_pt,
                             const unsigned& ieqn_local);
 
    void get_residuals(GeneralisedElement* const& elem_pt,
                       Vector<double>& residuals);
 
    void get_jacobian(GeneralisedElement* const& elem_pt,
                      Vector<double>& residuals,
                      DenseMatrix<double>& jacobian);
 
    void get_all_vectors_and_matrices(GeneralisedElement* const& elem_pt,
                                      Vector<Vector<double>>& vec,
                                      Vector<DenseMatrix<double>>& matrix);
 
    ~EigenProblemHandler() {}
  };
 
 
  //=============================================================
  //===============================================================
  class ParallelResidualsHandler : public AssemblyHandler
  {
    AssemblyHandler* Assembly_handler_pt;
 
  public:
    ParallelResidualsHandler(AssemblyHandler* const& assembly_handler_pt)
      : Assembly_handler_pt(assembly_handler_pt)
    {
    }
 
    unsigned ndof(GeneralisedElement* const& elem_pt)
    {
      return Assembly_handler_pt->ndof(elem_pt);
    }
 
    unsigned long eqn_number(GeneralisedElement* const& elem_pt,
                             const unsigned& ieqn_local)
    {
      return Assembly_handler_pt->eqn_number(elem_pt, ieqn_local);
    }
 
    void get_residuals(GeneralisedElement* const& elem_pt,
                       Vector<double>& residuals)
    {
      Assembly_handler_pt->get_residuals(elem_pt, residuals);
    }
 
 
    void get_jacobian(GeneralisedElement* const& elem_pt,
                      Vector<double>& residuals,
                      DenseMatrix<double>& jacobian)
    {
      Assembly_handler_pt->get_jacobian(elem_pt, residuals, jacobian);
    }
 
 
    void get_all_vectors_and_matrices(GeneralisedElement* const& elem_pt,
                                      Vector<Vector<double>>& vec,
                                      Vector<DenseMatrix<double>>& matrix)
    {
      Assembly_handler_pt->get_residuals(elem_pt, vec[0]);
    }
 
    ~ParallelResidualsHandler() {}
  };
 
 
  //==========================================================================
  //==========================================================================
  class ParameterDerivativeHandler : public AssemblyHandler
  {
    double* Parameter_pt;
 
    AssemblyHandler* Assembly_handler_pt;
 
  public:
    ParameterDerivativeHandler(AssemblyHandler* const& assembly_handler_pt,
                               double* const& parameter_pt)
      : Parameter_pt(parameter_pt), Assembly_handler_pt(assembly_handler_pt)
    {
    }
 
    unsigned ndof(GeneralisedElement* const& elem_pt)
    {
      return Assembly_handler_pt->ndof(elem_pt);
    }
 
    unsigned long eqn_number(GeneralisedElement* const& elem_pt,
                             const unsigned& ieqn_local)
    {
      return Assembly_handler_pt->eqn_number(elem_pt, ieqn_local);
    }
 
    void get_residuals(GeneralisedElement* const& elem_pt,
                       Vector<double>& residuals)
    {
      Assembly_handler_pt->get_dresiduals_dparameter(
        elem_pt, Parameter_pt, residuals);
    }
 
 
    void get_jacobian(GeneralisedElement* const& elem_pt,
                      Vector<double>& residuals,
                      DenseMatrix<double>& jacobian)
    {
      Assembly_handler_pt->get_djacobian_dparameter(
        elem_pt, Parameter_pt, residuals, jacobian);
    }
  };
 
 
  //========================================================================
  //========================================================================
  class AugmentedBlockFoldLinearSolver : public LinearSolver
  {
    LinearSolver* Linear_solver_pt;
 
    // Pointer to the problem, used in the resolve
    Problem* Problem_pt;
 
    DoubleVector* Alpha_pt;
 
    DoubleVector* E_pt;
 
  public:
    AugmentedBlockFoldLinearSolver(LinearSolver* const linear_solver_pt)
      : Linear_solver_pt(linear_solver_pt), Problem_pt(0), Alpha_pt(0), E_pt(0)
    {
    }
 
    ~AugmentedBlockFoldLinearSolver();
 
    void solve(Problem* const& problem_pt, DoubleVector& result);
 
    void solve(DoubleMatrixBase* const& matrix_pt,
               const DoubleVector& rhs,
               DoubleVector& result)
    {
      throw OomphLibError(
        "Linear-algebra interface does not make sense for this linear solver\n",
        OOMPH_CURRENT_FUNCTION,
        OOMPH_EXCEPTION_LOCATION);
    }
 
    void solve(DoubleMatrixBase* const& matrix_pt,
               const Vector<double>& rhs,
               Vector<double>& result)
    {
      throw OomphLibError(
        "Linear-algebra interface does not make sense for this linear solver\n",
        OOMPH_CURRENT_FUNCTION,
        OOMPH_EXCEPTION_LOCATION);
    }
 
    void resolve(const DoubleVector& rhs, DoubleVector& result);
 
    LinearSolver* linear_solver_pt() const
    {
      return Linear_solver_pt;
    }
  };
 
 
  //===================================================================
  //======================================================================
  class FoldHandler : public AssemblyHandler
  {
    friend class AugmentedBlockFoldLinearSolver;
 
    enum
    {
      Full_augmented,
      Block_J,
      Block_augmented_J
    };
 
    unsigned Solve_which_system;
 
    Problem* Problem_pt;
 
    unsigned Ndof;
 
    Vector<double> Phi;
 
    Vector<double> Y;
 
    Vector<int> Count;
 
    double* Parameter_pt;
 
  public:
    FoldHandler(Problem* const& problem_pt, double* const& parameter_pt);
 
 
    FoldHandler(Problem* const& problem_pt,
                double* const& parameter_pt,
                const DoubleVector& eigenvector);
 
    FoldHandler(Problem* const& problem_pt,
                double* const& parameter_pt,
                const DoubleVector& eigenvector,
                const DoubleVector& normalisation);
 
 
    ~FoldHandler();
 
    unsigned ndof(GeneralisedElement* const& elem_pt);
 
    unsigned long eqn_number(GeneralisedElement* const& elem_pt,
                             const unsigned& ieqn_local);
 
    void get_residuals(GeneralisedElement* const& elem_pt,
                       Vector<double>& residuals);
 
    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
    {
      return 1;
    }
 
    double* bifurcation_parameter_pt() const
    {
      return Parameter_pt;
    }
 
    void get_eigenfunction(Vector<DoubleVector>& eigenfunction);
 
    void solve_augmented_block_system();
 
    void solve_block_system();
 
    void solve_full_system();
  };
 
 
  //========================================================================
  //========================================================================
  class BlockPitchForkLinearSolver : public LinearSolver
  {
    LinearSolver* Linear_solver_pt;
 
    Problem* Problem_pt;
 
    DoubleVector* B_pt;
 
    DoubleVector* C_pt;
 
    DoubleVector* D_pt;
 
    DoubleVector* dJy_dparam_pt;
 
  public:
    BlockPitchForkLinearSolver(LinearSolver* const linear_solver_pt)
      : Linear_solver_pt(linear_solver_pt),
        Problem_pt(0),
        B_pt(0),
        C_pt(0),
        D_pt(0),
        dJy_dparam_pt(0)
    {
    }
 
    ~BlockPitchForkLinearSolver();
 
    void solve(Problem* const& problem_pt, DoubleVector& result);
 
    void solve(DoubleMatrixBase* const& matrix_pt,
               const DoubleVector& rhs,
               DoubleVector& result)
    {
      throw OomphLibError(
        "Linear-algebra interface does not make sense for this linear solver\n",
        OOMPH_CURRENT_FUNCTION,
        OOMPH_EXCEPTION_LOCATION);
    }
 
    void solve(DoubleMatrixBase* const& matrix_pt,
               const Vector<double>& rhs,
               Vector<double>& result)
    {
      throw OomphLibError(
        "Linear-algebra interface does not make sense for this linear solver\n",
        OOMPH_CURRENT_FUNCTION,
        OOMPH_EXCEPTION_LOCATION);
    }
 
 
    void resolve(const DoubleVector& rhs, DoubleVector& result);
 
    LinearSolver* linear_solver_pt() const
    {
      return Linear_solver_pt;
    }
  };
 
 
  //========================================================================
  //========================================================================
  class AugmentedBlockPitchForkLinearSolver : public LinearSolver
  {
    LinearSolver* Linear_solver_pt;
 
    Problem* Problem_pt;
 
    DoubleVector* Alpha_pt;
 
    DoubleVector* E_pt;
 
  public:
    AugmentedBlockPitchForkLinearSolver(LinearSolver* const linear_solver_pt)
      : Linear_solver_pt(linear_solver_pt), Problem_pt(0), Alpha_pt(0), E_pt(0)
    {
    }
 
    ~AugmentedBlockPitchForkLinearSolver();
 
    void solve(Problem* const& problem_pt, DoubleVector& result);
 
    void solve(DoubleMatrixBase* const& matrix_pt,
               const DoubleVector& rhs,
               DoubleVector& result)
    {
      throw OomphLibError(
        "Linear-algebra interface does not make sense for this linear solver\n",
        OOMPH_CURRENT_FUNCTION,
        OOMPH_EXCEPTION_LOCATION);
    }
 
    void solve(DoubleMatrixBase* const& matrix_pt,
               const Vector<double>& rhs,
               Vector<double>& result)
    {
      throw OomphLibError(
        "Linear-algebra interface does not make sense for this linear solver\n",
        OOMPH_CURRENT_FUNCTION,
        OOMPH_EXCEPTION_LOCATION);
    }
 
 
    void resolve(const DoubleVector& rhs, DoubleVector& result);
 
    LinearSolver* linear_solver_pt() const
    {
      return Linear_solver_pt;
    }
  };
 
  //========================================================================
  //========================================================================
  class PitchForkHandler : public AssemblyHandler
  {
    friend class BlockPitchForkLinearSolver;
    friend class AugmentedBlockPitchForkLinearSolver;
 
    enum
    {
      Full_augmented,
      Block_J,
      Block_augmented_J
    };
 
    unsigned Solve_which_system;
 
    Problem* Problem_pt;
 
    AssemblyHandler* Assembly_handler_pt;
 
    unsigned Ndof;
 
    LinearAlgebraDistribution* Dof_distribution_pt;
 
    LinearAlgebraDistribution* Augmented_dof_distribution_pt;
 
    double Sigma;
 
    DoubleVectorWithHaloEntries Y;
 
    DoubleVectorWithHaloEntries Psi;
 
    DoubleVectorWithHaloEntries C;
 
    DoubleVectorWithHaloEntries Count;
 
    Vector<unsigned> Global_eqn_number;
 
    double* Parameter_pt;
 
    // The total number of elements in the problem
    unsigned Nelement;
 
#ifdef OOMPH_HAS_MPI
    bool Distributed;
#endif
 
    inline unsigned global_eqn_number(const unsigned& i)
    {
#ifdef OOMPH_HAS_MPI
      // If the problem is distributed I have to do something
      if (Distributed)
      {
        return Global_eqn_number[i];
      }
      // Otherwise it's just i
      else
#endif
      {
        return i;
      }
    }
 
  public:
    PitchForkHandler(Problem* const& problem_pt,
                     AssemblyHandler* const& assembly_handler_pt,
                     double* const& parameter_pt,
                     const DoubleVector& symmetry_vector);
 
    ~PitchForkHandler();
 
    // Has this been called
 
    unsigned ndof(GeneralisedElement* const& elem_pt);
 
    unsigned long eqn_number(GeneralisedElement* const& elem_pt,
                             const unsigned& ieqn_local);
 
    void get_residuals(GeneralisedElement* const& elem_pt,
                       Vector<double>& residuals);
 
    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
    {
      return 2;
    }
 
    double* bifurcation_parameter_pt() const
    {
      return Parameter_pt;
    }
 
    void get_eigenfunction(Vector<DoubleVector>& eigenfunction);
 
#ifdef OOMPH_HAS_MPI
    void synchronise();
#endif
 
 
    void solve_augmented_block_system();
 
    void solve_block_system();
 
    void solve_full_system();
  };
 
  //========================================================================
  //========================================================================
  class BlockHopfLinearSolver : public LinearSolver
  {
    LinearSolver* Linear_solver_pt;
 
    Problem* Problem_pt;
 
    DoubleVector* A_pt;
 
    DoubleVector* E_pt;
 
    DoubleVector* G_pt;
 
  public:
    BlockHopfLinearSolver(LinearSolver* const linear_solver_pt)
      : Linear_solver_pt(linear_solver_pt),
        Problem_pt(0),
        A_pt(0),
        E_pt(0),
        G_pt(0)
    {
    }
 
    ~BlockHopfLinearSolver();
 
    void solve_for_two_rhs(Problem* const& problem_pt,
                           DoubleVector& result,
                           const DoubleVector& rhs2,
                           DoubleVector& result2);
 
    void solve(Problem* const& problem_pt, DoubleVector& result);
 
    void solve(DoubleMatrixBase* const& matrix_pt,
               const DoubleVector& rhs,
               DoubleVector& result)
    {
      throw OomphLibError(
        "Linear-algebra interface does not make sense for this linear solver\n",
        OOMPH_CURRENT_FUNCTION,
        OOMPH_EXCEPTION_LOCATION);
    }
 
    void solve(DoubleMatrixBase* const& matrix_pt,
               const Vector<double>& rhs,
               Vector<double>& result)
    {
      throw OomphLibError(
        "Linear-algebra interface does not make sense for this linear solver\n",
        OOMPH_CURRENT_FUNCTION,
        OOMPH_EXCEPTION_LOCATION);
    }
 
 
    void resolve(const DoubleVector& rhs, DoubleVector& result);
 
    LinearSolver* linear_solver_pt() const
    {
      return Linear_solver_pt;
    }
  };
 
 
  //===============================================================
  //===========================================================================
  class HopfHandler : public AssemblyHandler
  {
    friend class BlockHopfLinearSolver;
 
    unsigned Solve_which_system;
 
    Problem* Problem_pt;
 
    double* Parameter_pt;
 
    unsigned Ndof;
 
    double Omega;
 
    Vector<double> Phi;
 
    Vector<double> Psi;
 
    Vector<double> C;
 
    Vector<int> Count;
 
  public:
    HopfHandler(Problem* const& problem_pt, double* const& parameter_pt);
 
    HopfHandler(Problem* const& problem_pt,
                double* const& paramter_pt,
                const double& omega,
                const DoubleVector& phi,
                const DoubleVector& psi);
 
    ~HopfHandler();
 
    unsigned ndof(GeneralisedElement* const& elem_pt);
 
    unsigned long eqn_number(GeneralisedElement* const& elem_pt,
                             const unsigned& ieqn_local);
 
    void get_residuals(GeneralisedElement* const& elem_pt,
                       Vector<double>& residuals);
 
    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
    {
      return 3;
    }
 
    double* bifurcation_parameter_pt() const
    {
      return Parameter_pt;
    }
 
    void get_eigenfunction(Vector<DoubleVector>& eigenfunction);
 
    const double& omega() const
    {
      return Omega;
    }
 
    void solve_standard_system();
 
    void solve_complex_system();
 
    void solve_full_system();
  };
 
 
} // namespace oomph
#endif