// 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-2026 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//==================================================================== #ifndef OOMPH_PSEUDO_ELASTIC_SUBSIDIARY_PRECONDITIONER #define OOMPH_PSEUDO_ELASTIC_SUBSIDIARY_PRECONDITIONER // includes #include "generic/problem.h" #include "generic/block_preconditioner.h" #include "generic/preconditioner.h" #include "generic/SuperLU_preconditioner.h" #include "generic/matrix_vector_product.h" #include "generic/general_purpose_preconditioners.h" #include "generic/general_purpose_block_preconditioners.h" #ifdef OOMPH_HAS_HYPRE #include "generic/hypre_solver.h" #endif #ifdef OOMPH_HAS_TRILINOS #include "generic/trilinos_solver.h" #endif namespace oomph { //============================================================================= /// Functions to create instances of optimal subsidiary operators for /// the PseudoElasticPreconditioner. By default we use hypre for the /// the elastic blocks but can use Trilinos ML too. //============================================================================= namespace Pseudo_Elastic_Preconditioner_Subsidiary_Operator_Helper { #ifdef OOMPH_HAS_HYPRE /// Hypre AMG w/ GS smoothing for the augmented elastic /// subsidiary linear systems Preconditioner* get_elastic_preconditioner_hypre(); /// AMG w/ GS smoothing for the augmented elastic subsidiary linear /// systems -- calls Hypre version to stay consistent with previous default Preconditioner* get_elastic_preconditioner(); #endif #ifdef OOMPH_HAS_TRILINOS /// TrilinosML smoothing for the augmented elastic /// subsidiary linear systems Preconditioner* get_elastic_preconditioner_trilinos_ml(); /// CG with diagonal preconditioner for the lagrange multiplier /// subsidiary linear systems. Preconditioner* get_lagrange_multiplier_preconditioner(); #endif } // namespace Pseudo_Elastic_Preconditioner_Subsidiary_Operator_Helper /////////////////////////////////////////////////////////////////////////////// /////////////////////////////////////////////////////////////////////////////// /////////////////////////////////////////////////////////////////////////////// //============================================================================= /// A subsidiary preconditioner for the pseudo-elastic FSI /// preconditioner. Also a stand-alone preconditioner for the problem of /// non-linear elasticity subject to prescribed displacement by Lagrange /// multiplier. /// \b Enumeration of Elastic DOF types in the Pseudo-Elastic Elements /// The method get_dof_types_for_unknowns() must be implemented such that /// DOFs subject be Lagrange multiplier and DOFs NOT subject to Lagrange /// multiplier have different labels. For example in a 3D problem there are /// 6 DOF types and the following labelling must be implemented: /// 0 - x displacement (without lagr mult traction) /// 1 - y displacement (without lagr mult traction) /// 2 - z displacement (without lagr mult traction) /// 4 - x displacement (with lagr mult traction) /// 5 - y displacement (with lagr mult traction) /// 6 - z displacement (with lagr mult traction) //============================================================================= class PseudoElasticPreconditioner : public BlockPreconditioner { /// PseudoElasticFSIPreconditioner is a friend to access the private /// *_preconditioner_solve(...) method friend class PseudoElasticFSIPreconditioner; public: /// This preconditioner includes the option to use subsidiary /// operators other than ExactPreconditioner for this problem. /// This is the typedef of a function that should return an instance /// of a subsidiary preconditioning operator. This preconditioner is /// responsible for the destruction of the subsidiary preconditioners. typedef Preconditioner* (*SubsidiaryPreconditionerFctPt)(); /// The augmented elasticity system can be preconditioned in one /// of four ways. /// 0 - Exact preconditioner /// 1 - Block diagonal preconditioning /// 2 - Block upper triangular preconditioner /// 3 - Block lower triangular preconditioner /// We group together the different components of the displacement vector /// field for the block decomposition. enum Elastic_preconditioner_type { Exact_block_preconditioner, Block_diagonal_preconditioner, Block_lower_triangular_preconditioner, Block_upper_triangular_preconditioner }; /// Default (and only) constructor. PseudoElasticPreconditioner() { // null pointers Lagrange_multiplier_subsidiary_preconditioner_function_pt = 0; Elastic_subsidiary_preconditioner_function_pt = 0; Elastic_preconditioner_pt = 0; // set defaults Use_inf_norm_of_s_scaling = true; E_preconditioner_type = Exact_block_preconditioner; // resize the Mesh_pt Lagrange_multiplier_mesh_pt = 0; Elastic_mesh_pt = 0; } /// destructor virtual ~PseudoElasticPreconditioner() { this->clean_up_memory(); } /// Broken copy constructor PseudoElasticPreconditioner(const PseudoElasticPreconditioner&) = delete; /// Broken assignment operator // Commented out broken assignment operator because this can lead to a // conflict warning when used in the virtual inheritence hierarchy. // Essentially the compiler doesn't realise that two separate // implementations of the broken function are the same and so, quite // rightly, it shouts. /*void operator= (const PseudoElasticPreconditioner&) = delete;*/ /// Setup method for the PseudoElasticPreconditioner. void setup(); /// Apply the preconditioner. Method implemented in two /// other methods (elastic and lagrange multiplier subsidiary /// preocnditioner) for the PseudoElasticFSIPreconditioner void preconditioner_solve(const DoubleVector& r, DoubleVector& z) { this->elastic_preconditioner_solve(r, z); this->lagrange_multiplier_preconditioner_solve(r, z); } /// Access function to mesh containing the block-preconditionable /// elastic elements void set_elastic_mesh(Mesh* mesh_pt) { Elastic_mesh_pt = mesh_pt; } /// Access function to mesh containing the block-preconditionable /// lagrange multiplier elements void set_lagrange_multiplier_mesh(Mesh* mesh_pt) { Lagrange_multiplier_mesh_pt = mesh_pt; } /// Call to use the inf norm of S as scaling void enable_inf_norm_of_s_scaling() { Use_inf_norm_of_s_scaling = true; } /// Call to use no scaling void disable_inf_norm_of_s_scaling() { Use_inf_norm_of_s_scaling = false; } /// By default the Lagrange multiplier subsidiary systems are /// preconditioner with ExactPreconditioner. For a different /// preconditioner, pass a function to this /// method returning a different subsidiary operator. void set_lagrange_multiplier_subsidiary_preconditioner( SubsidiaryPreconditionerFctPt prec_fn) { Lagrange_multiplier_subsidiary_preconditioner_function_pt = prec_fn; } /// By default the elastic subsidiary systems are /// preconditioner with ExactPreconditioner. For a different /// preconditioner, pass a function to this /// method returning a different subsidiary operator. void set_elastic_subsidiary_preconditioner( SubsidiaryPreconditionerFctPt prec_fn) { Elastic_subsidiary_preconditioner_function_pt = prec_fn; } /// Set the type of preconditioner applied to the elastic: /// 0 - Exact preconditioner /// 1 - Block diagonal preconditioning /// 2 - Block upper triangular preconditioner /// 3 - Block lower triangular preconditioner /// We group together the different components of the displacement vector /// field for the block decomposition. Elastic_preconditioner_type& elastic_preconditioner_type() { return E_preconditioner_type; } /// Clears the memory. void clean_up_memory(); private: /// Apply the elastic subsidiary preconditioner. void elastic_preconditioner_solve(const DoubleVector& r, DoubleVector& z); /// Apply the lagrange multiplier subsidiary preconditioner. void lagrange_multiplier_preconditioner_solve(const DoubleVector& r, DoubleVector& z); /// The scaling. Defaults to infinity norm of S. double Scaling; /// boolean indicating whether the inf-norm of S should be used as /// scaling. Default = true; bool Use_inf_norm_of_s_scaling; /// An unsigned indicating which method should be used for /// preconditioning the solid component. Elastic_preconditioner_type E_preconditioner_type; /// the dimension of the problem unsigned Dim; /// storage for the preconditioner for the solid system Preconditioner* Elastic_preconditioner_pt; /// lagrange multiplier preconditioner pt Vector Lagrange_multiplier_preconditioner_pt; /// The Lagrange multiplier subsidiary preconditioner function pointer SubsidiaryPreconditionerFctPt Lagrange_multiplier_subsidiary_preconditioner_function_pt; /// The solid subsidiary preconditioner function pointer SubsidiaryPreconditionerFctPt Elastic_subsidiary_preconditioner_function_pt; /// Pointer to the mesh containing the solid elements Mesh* Elastic_mesh_pt; /// Pointer to the mesh containing the Lagrange multiplier elements Mesh* Lagrange_multiplier_mesh_pt; }; /////////////////////////////////////////////////////////////////////////////// /////////////////////////////////////////////////////////////////////////////// /////////////////////////////////////////////////////////////////////////////// //============================================================================= /// A subsidiary preconditioner for the pseudo-elastic FSI /// preconditioner. Also a stand-alone preconditioner for the problem of /// non-linear elasticity subject to prescribed displacement by Lagrange /// multiplier.. /// \b Enumeration of Elastic DOF types in the Pseudo-Elastic Elements /// The method get_dof_types_for_unknowns() must be implemented such that /// DOFs subject be Lagrange multiplier and DOFs NOT subject to Lagrange /// multiplier have different labels. For example in a 3D problem there are /// 6 DOF types and the following labelling must be implemented: /// 0 - x displacement (without lagr mult traction) /// 1 - y displacement (without lagr mult traction) /// 2 - z displacement (without lagr mult traction) /// 4 - x displacement (with lagr mult traction) /// 5 - y displacement (with lagr mult traction) /// 6 - z displacement (with lagr mult traction) //============================================================================= class PseudoElasticPreconditionerOld : public BlockPreconditioner { /// PseudoElasticFSIPreconditioner is a friend to access the private /// *_preconditioner_solve(...) method friend class PseudoElasticFSIPreconditioner; public: /// This preconditioner includes the option to use subsidiary /// operators other than ExactPreconditioner for this problem. /// This is the typedef of a function that should return an instance /// of a subsidiary preconditioning operator. This preconditioner is /// responsible for the destruction of the subsidiary preconditioners. typedef Preconditioner* (*SubsidiaryPreconditionerFctPt)(); /// The augmented elasticity system can be preconditioned in one /// of four ways. /// 0 - Exact preconditioner /// 1 - Block diagonal preconditioning /// 2 - Block upper triangular preconditioner /// 3 - Block lower triangular preconditioner /// We group together the different components of the displacement vector /// field for the block decomposition. enum Elastic_preconditioner_type { Exact_block_preconditioner, Block_diagonal_preconditioner, Block_lower_triangular_preconditioner, Block_upper_triangular_preconditioner }; /// Default (and only) constructor. PseudoElasticPreconditionerOld() { // null pointers Lagrange_multiplier_subsidiary_preconditioner_function_pt = 0; Elastic_subsidiary_preconditioner_function_pt = 0; Elastic_preconditioner_pt = 0; // set defaults Use_inf_norm_of_s_scaling = true; E_preconditioner_type = Exact_block_preconditioner; // resize the Mesh_pt Lagrange_multiplier_mesh_pt = 0; Elastic_mesh_pt = 0; } /// destructor virtual ~PseudoElasticPreconditionerOld() { this->clean_up_memory(); } /// Broken copy constructor PseudoElasticPreconditionerOld(const PseudoElasticPreconditionerOld&) = delete; /// Broken assignment operator /*void operator=(const PseudoElasticPreconditionerOld&) = delete;*/ /// Setup method for the PseudoElasticPreconditionerOld. void setup(); /// Apply the preconditioner. Method implemented in two /// other methods (elastic and lagrange multiplier subsidiary /// preocnditioner) for the PseudoElasticFSIPreconditioner void preconditioner_solve(const DoubleVector& r, DoubleVector& z) { this->elastic_preconditioner_solve(r, z); this->lagrange_multiplier_preconditioner_solve(r, z); } /// Access function to mesh containing the block-preconditionable /// elastic elements void set_elastic_mesh(Mesh* mesh_pt) { Elastic_mesh_pt = mesh_pt; } /// Access function to mesh containing the block-preconditionable /// lagrange multiplier elements void set_lagrange_multiplier_mesh(Mesh* mesh_pt) { Lagrange_multiplier_mesh_pt = mesh_pt; } /// Call to use the inf norm of S as scaling void enable_inf_norm_of_s_scaling() { Use_inf_norm_of_s_scaling = true; } /// Call to use no scaling void disable_inf_norm_of_s_scaling() { Use_inf_norm_of_s_scaling = false; } /// By default the Lagrange multiplier subsidiary systems are /// preconditioner with ExactPreconditioner. For a different /// preconditioner, pass a function to this /// method returning a different subsidiary operator. void set_lagrange_multiplier_subsidiary_preconditioner( SubsidiaryPreconditionerFctPt prec_fn) { Lagrange_multiplier_subsidiary_preconditioner_function_pt = prec_fn; } /// By default the elastic subsidiary systems are /// preconditioner with ExactPreconditioner. For a different /// preconditioner, pass a function to this /// method returning a different subsidiary operator. void set_elastic_subsidiary_preconditioner( SubsidiaryPreconditionerFctPt prec_fn) { Elastic_subsidiary_preconditioner_function_pt = prec_fn; } /// Set the type of preconditioner applied to the elastic: /// 0 - Exact preconditioner /// 1 - Block diagonal preconditioning /// 2 - Block upper triangular preconditioner /// 3 - Block lower triangular preconditioner /// We group together the different components of the displacement vector /// field for the block decomposition. Elastic_preconditioner_type& elastic_preconditioner_type() { return E_preconditioner_type; } /// Clears the memory. void clean_up_memory(); private: /// Apply the elastic subsidiary preconditioner. void elastic_preconditioner_solve(const DoubleVector& r, DoubleVector& z); /// Apply the lagrange multiplier subsidiary preconditioner. void lagrange_multiplier_preconditioner_solve(const DoubleVector& r, DoubleVector& z); /// The scaling. Defaults to infinity norm of S. double Scaling; /// boolean indicating whether the inf-norm of S should be used as /// scaling. Default = true; bool Use_inf_norm_of_s_scaling; /// An unsigned indicating which method should be used for /// preconditioning the solid component. Elastic_preconditioner_type E_preconditioner_type; /// the dimension of the problem unsigned Dim; /// storage for the preconditioner for the solid system Preconditioner* Elastic_preconditioner_pt; /// lagrange multiplier preconditioner pt Vector Lagrange_multiplier_preconditioner_pt; /// The Lagrange multiplier subsidary preconditioner function pointer SubsidiaryPreconditionerFctPt Lagrange_multiplier_subsidiary_preconditioner_function_pt; /// The solid subsidiary preconditioner function pointer SubsidiaryPreconditionerFctPt Elastic_subsidiary_preconditioner_function_pt; /// Pointer to the mesh containing the solid elements Mesh* Elastic_mesh_pt; /// Pointer to the mesh containing the Lagrange multiplier elements Mesh* Lagrange_multiplier_mesh_pt; }; //////////////////////////////////////////////////////////////////////////////// //////////////////////////////////////////////////////////////////////////////// //////////////////////////////////////////////////////////////////////////////// //============================================================================= /// Subsidiary helper preconditioner for the PseudoElasticPreconditioner. /// Required to construct the augmented elastic system prior to /// preconditioning. /// NOTE: /// 1. This is only intended to be used as a subsidiary preconditioner within /// the PseudoElasticPreconditioner. /// 2. If this preconditioner has N DOF types then the first N/2 are assumed /// to be ordinary solid DOF types, and the second N/2 are the solid DOF types /// with lagrange multiplier tractions applied. /// 3. By default this preconditioner uses a superlu preconditioner. //============================================================================= class PseudoElasticPreconditionerSubsidiaryPreconditionerOld : public BlockPreconditioner { public: /// typedef for a function that allows other preconditioners to be /// emplyed to solve the subsidiary linear systems. /// The function should return a pointer to the requred subsidiary /// preconditioner generated using new. This preconditioner is responsible /// for the destruction of the subsidiary preconditioners. typedef Preconditioner* (*SubsidiaryPreconditionerFctPt)(); /// Constructor PseudoElasticPreconditionerSubsidiaryPreconditionerOld() { Scaling = 1.0; Preconditioner_pt = 0; Subsidiary_preconditioner_function_pt = 0; } /// Destructor ~PseudoElasticPreconditionerSubsidiaryPreconditionerOld() { this->clean_up_memory(); } /// Broken copy constructor PseudoElasticPreconditionerSubsidiaryPreconditionerOld( const PseudoElasticPreconditionerSubsidiaryPreconditionerOld&) = delete; /// Broken assignment operator /*void operator=(const PseudoElasticPreconditionerSubsidiaryPreconditionerOld&) = delete;*/ // Setup the preconditioner void setup(); // Apply the preconditioner void preconditioner_solve(const DoubleVector& r, DoubleVector& z); /// Specify the scaling. Default is 1.0 Must be called before /// setup(...). double& scaling() { return Scaling; } /// access function to set the subsidiary preconditioner function. void set_subsidiary_preconditioner_function( SubsidiaryPreconditionerFctPt sub_prec_fn) { Subsidiary_preconditioner_function_pt = sub_prec_fn; } private: /// clears the memory void clean_up_memory() { delete Preconditioner_pt; Preconditioner_pt = 0; } // the augmentation scaling double Scaling; /// the preconditioner pt Preconditioner* Preconditioner_pt; /// the SubisidaryPreconditionerFctPt SubsidiaryPreconditionerFctPt Subsidiary_preconditioner_function_pt; }; // end of PseudoElasticPreconditionerSubsidiaryPreconditionerOld //////////////////////////////////////////////////////////////////////////////// //////////////////////////////////////////////////////////////////////////////// //////////////////////////////////////////////////////////////////////////////// //============================================================================= /// Subsidiary helper preconditioner for the PseudoElasticPreconditioner. /// Required for block preconditioner of the augmented elastic subsidiary /// problem. /// NOTE: /// 1. This is only intended to be used as a subsidiary preconditioner within /// the PseudoElasticPreconditioner. /// 2. If this preconditioner has N DOF types then the first N/2 are assumed /// to be ordinary solid DOF types, and the second N/2 are the solid DOF types /// with lagrange multiplier tractions applied. /// 3. By default this preconditioner uses a superlu preconditioner. //============================================================================= class PseudoElasticPreconditionerSubsidiaryBlockPreconditionerOld : public BlockPreconditioner { public: /// This preconditioner includes the option to use subsidiary /// operators other than ExactPreconditioner for this problem. /// This is the typedef of a function that should return an instance /// of a subsidiary preconditioning operator. This preconditioner is /// responsible for the destruction of the subsidiary preconditioners. typedef Preconditioner* (*SubsidiaryPreconditionerFctPt)(); /// Constructor. (By default this preconditioner is upper triangular). PseudoElasticPreconditionerSubsidiaryBlockPreconditionerOld() : BlockPreconditioner() { // null the subsidiary preconditioner function pointer Subsidiary_preconditioner_function_pt = 0; // default to block diagonal Method = 0; // default scaling = 1.0 Scaling = 1.0; }; /// Destructor ~PseudoElasticPreconditionerSubsidiaryBlockPreconditionerOld() { this->clean_up_memory(); } /// Broken copy constructor PseudoElasticPreconditionerSubsidiaryBlockPreconditionerOld( const PseudoElasticPreconditionerSubsidiaryBlockPreconditionerOld&) = delete; /// Broken assignment operator /*void operator= (const PseudoElasticPreconditionerSubsidiaryBlockPreconditionerOld&) = delete;*/ /// clean up the memory void clean_up_memory(); /// Setup the preconditioner void setup(); /// Apply preconditioner to r void preconditioner_solve(const DoubleVector& res, DoubleVector& z); /// access function to set the subsidiary preconditioner function. void set_subsidiary_preconditioner_function( SubsidiaryPreconditionerFctPt sub_prec_fn) { Subsidiary_preconditioner_function_pt = sub_prec_fn; }; /// use as a block diagonal preconditioner void use_block_diagonal_approximation() { Method = 0; } /// Use as an upper triangular preconditioner void use_upper_triangular_approximation() { Method = 1; } /// Use as a lower triangular preconditioner void use_lower_triangular_approximation() { Method = 2; } /// Specify the scaling. Default is 1.0 Must be set before /// setup(...). double& scaling() { return Scaling; } private: /// Vector of SuperLU preconditioner pointers for storing the /// preconditioners for each diagonal block Vector Diagonal_block_preconditioner_pt; /// Matrix of matrix vector product operators for the off diagonals DenseMatrix Off_diagonal_matrix_vector_products; /// the preconditioning method. /// 0 - block diagonal /// 1 - upper triangular /// 2 - lower triangular unsigned Method; /// The SubisidaryPreconditionerFctPt SubsidiaryPreconditionerFctPt Subsidiary_preconditioner_function_pt; /// The scaling. default 1.0. double Scaling; }; /////////////////////////////////////////////////////////////////////////////// /////////////////////////////////////////////////////////////////////////////// /////////////////////////////////////////////////////////////////////////////// // /* //============================================================================= /// A helper class for PseudoElasticPreconditioner. /// Note that this is NOT actually a functioning preconditioner. /// We simply derive from this class to get access to the blocks. //============================================================================= class PseudoElasticPreconditionerScalingHelperOld : public BlockPreconditioner { public: /// The constructor. /// NOTE: /// 1. master_prec_pt should point to the /// PseudoElasticPreconditioner. /// 2. matrix_pt should point to the jacobian. /// 3. The vector dof_list should contain the full list of /// DOFS associated with the solid subsidiary system. /// 4. "solid_mesh_pt" should be a pointer to the solid mesh used in the /// master preconditioner. PseudoElasticPreconditionerScalingHelperOld( BlockPreconditioner* master_prec_pt, CRDoubleMatrix* matrix_pt, Vector& dof_list, const Mesh* const solid_mesh_pt, const OomphCommunicator* comm_pt) { // turn into a subisiary preconditioner this->turn_into_subsidiary_block_preconditioner(master_prec_pt, dof_list); // all dofs are of the same block type Vector dof_to_block_map(dof_list.size(), 0); // store the matrix_pt set_matrix_pt(matrix_pt); // set the mesh this->set_nmesh(1); this->set_mesh(0, solid_mesh_pt); // set the communicator pointer this->set_comm_pt(comm_pt); // call block_setup(...) this->block_setup(dof_to_block_map); } /// Destructor. ~PseudoElasticPreconditionerScalingHelperOld() { this->clear_block_preconditioner_base(); } /// Broken copy constructor PseudoElasticPreconditionerScalingHelperOld( const PseudoElasticPreconditionerScalingHelperOld&) = delete; /// Broken assignment operator /*void operator=(const PseudoElasticPreconditionerScalingHelperOld&) = * delete;*/ /// returns the infinite norm of S double s_inf_norm() { CRDoubleMatrix* m_pt = new CRDoubleMatrix; this->get_block(0, 0, *m_pt); double s_inf_norm = m_pt->inf_norm(); delete m_pt; return s_inf_norm; } // broken preconditioner setup void setup() { std::ostringstream error_message; error_message << "This method is intentionally broken. This class is not " "a functioning " << "preconditioner."; throw OomphLibError( error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION); } // broken preconditioner solve void preconditioner_solve(const DoubleVector& r, DoubleVector& z) { std::ostringstream error_message; error_message << "This method is intentionally broken. This class is not " "a functioning " << "preconditioner."; throw OomphLibError( error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION); } }; // end of PseudoElasticPreconditionerScalingHelperOld // */ } // namespace oomph #endif