ChasteGuides/ConvertingAssemblersToUseNewDesign

Converting assemblers to use the new FEM solver design

This page illustrates the simplest way to convert an assembler which uses the old assembler hierarchy to use the new design. See FiniteElementAssemblersAndSolvers for more information about the new design.

Below is the old SimpleLinearEllipticAssembler (a solver of linear elliptic PDEs using the old design), with all comments removed and new comments added to explain what would need to be changed in order to convert it to use the new design. Following this, the new solver SimpleLinearEllipticSolver which uses the new design is given.

The code below is also attached to the page.

Old SimpleLinearEllipticAssembler with comments on how to convert it 

#include "AbstractLinearEllipticPde.hpp"
// the following 3 headers are no longer needed
#include "AbstractLinearAssembler.hpp"
#include <boost/mpl/void.hpp>
#include <boost/mpl/if.hpp>

//  The CONCRETE template parameter can be deleted 
template<unsigned ELEMENT_DIM, unsigned SPACE_DIM, class CONCRETE = boost::mpl::void_> 
class SimpleLinearEllipticAssembler  
  : public AbstractLinearAssembler<ELEMENT_DIM, SPACE_DIM, 1, true, SimpleLinearEllipticAssembler<ELEMENT_DIM, SPACE_DIM, CONCRETE> >
      //  The concrete class will now inherit from two classes: AbstractAssemblerSolverHybrid<ELEMENT_DIM,SPACE_DIM,your_problem_dim,NORMAL>
      //    AND ONE OF:
      //      AbstractStaticLinearPdeSolver<ELEMENT_DIM,SPACE_DIM,your_problem_dim> 
      //    OR
      //      AbstractDynamicLinearPdeSolver<ELEMENT_DIM,SPACE_DIM,your_problem_dim> 
{
public:
    // All these (the following 6 lines) will typically no longer be needed 
    // and can be deleted
    static const unsigned E_DIM = ELEMENT_DIM; 
    static const unsigned S_DIM = SPACE_DIM; 
    static const unsigned P_DIM = 1u; 
    typedef SimpleLinearEllipticAssembler<ELEMENT_DIM, SPACE_DIM, CONCRETE> SelfType;
    typedef AbstractLinearAssembler<ELEMENT_DIM, SPACE_DIM, 1, true, SelfType> BaseClassType; 
    friend class AbstractStaticAssembler<ELEMENT_DIM, SPACE_DIM, 1, true, SelfType>;

protected:

    AbstractLinearEllipticPde<ELEMENT_DIM,SPACE_DIM>* mpEllipticPde;


    // ComputeMatrixTerm, ComputeVectorTerm and ComputeVectorSurfaceTerm don't need to be changed
    virtual c_matrix<double, 1*(ELEMENT_DIM+1), 1*(ELEMENT_DIM+1)> ComputeMatrixTerm(
        c_vector<double, ELEMENT_DIM+1>& rPhi,
        c_matrix<double, SPACE_DIM, ELEMENT_DIM+1>& rGradPhi,
        ChastePoint<SPACE_DIM>& rX,
        c_vector<double,1>& rU,
        c_matrix<double,1,SPACE_DIM>& rGradU,
        Element<ELEMENT_DIM,SPACE_DIM>* pElement)
    {
        ..
    }


    // ComputeMatrixTerm, ComputeVectorTerm and ComputeVectorSurfaceTerm don't need to be changed
    virtual c_vector<double,1*(ELEMENT_DIM+1)> ComputeVectorTerm(
        c_vector<double, ELEMENT_DIM+1>& rPhi,
        c_matrix<double, SPACE_DIM, ELEMENT_DIM+1>& rGradPhi,
        ChastePoint<SPACE_DIM>& rX,
        c_vector<double,1>& rU,
        c_matrix<double,1,SPACE_DIM>& rGradU,
        Element<ELEMENT_DIM,SPACE_DIM>* pElement)
    {
        ..
    }


    // ComputeMatrixTerm, ComputeVectorTerm and ComputeVectorSurfaceTerm don't need to be changed
    virtual c_vector<double, ELEMENT_DIM> ComputeVectorSurfaceTerm(const BoundaryElement<ELEMENT_DIM-1,SPACE_DIM>& rSurfaceElement,
            c_vector<double, ELEMENT_DIM>& rPhi,
            ChastePoint<SPACE_DIM>& rX)
    {
        ..
    }


public:
    // Constructor needs to call base class constructors of new parents (see new version below)
    SimpleLinearEllipticAssembler(AbstractTetrahedralMesh<ELEMENT_DIM,SPACE_DIM>* pMesh,
                                  AbstractLinearEllipticPde<ELEMENT_DIM,SPACE_DIM>* pPde,
                                  BoundaryConditionsContainer<ELEMENT_DIM,SPACE_DIM,1>* pBoundaryConditions,
                                  unsigned numQuadPoints = 2)
        : AbstractAssembler<ELEMENT_DIM,SPACE_DIM,1>(),
          BaseClassType(numQuadPoints)
    {
        mpEllipticPde = pPde;

        // The following two lines aren't valid anymore and should be deleted
        this->SetMesh(pMesh); 
        this->SetBoundaryConditionsContainer(pBoundaryConditions);
    }

    // PrepareForSolve() doesn't exist anymore. Use InitialiseForSolve if there's anything that needs to be done before the solve
    void PrepareForSolve()
    {
        BaseClassType::PrepareForSolve();
        assert(mpEllipticPde != NULL);
    }
    
    // NOTE: if you have PrepareForAssembleSystem or FinaliseAssembleSystem
    // implemented, these need to be change to PrepareForSetupLinearSystem(Vec currentSolution)
    // and  FinaliseLinearSystem(Vec currentSolution). 
    // 
    // The current time is not passed into either of these. If you need to use it do:
    // double time = PdeSimulationTime::GetTime()


    void InitialiseForSolve(Vec initialSolution)
    {
        // BaseClassType should be given as AbstractLinearPdeSolver
        BaseClassType::InitialiseForSolve(initialSolution);
        assert(this->mpLinearSystem);
        this->mpLinearSystem->SetMatrixIsSymmetric(true);
        this->mpLinearSystem->SetKspType("cg");
    }
};


// All this stuff can be deleted
template<unsigned ELEMENT_DIM, unsigned SPACE_DIM, class CONCRETE>
struct AssemblerTraits<SimpleLinearEllipticAssembler<ELEMENT_DIM, SPACE_DIM, CONCRETE> >
{
    typedef typename boost::mpl::if_<boost::mpl::is_void_<CONCRETE>,
                                     SimpleLinearEllipticAssembler<ELEMENT_DIM, SPACE_DIM, CONCRETE>,
                                     typename AssemblerTraits<CONCRETE>::CVT_CLASS>::type
            CVT_CLASS;
    typedef typename boost::mpl::if_<boost::mpl::is_void_<CONCRETE>,
                                     SimpleLinearEllipticAssembler<ELEMENT_DIM, SPACE_DIM, CONCRETE>,
                                     typename AssemblerTraits<CONCRETE>::CMT_CLASS>::type
            CMT_CLASS;
    typedef typename boost::mpl::if_<boost::mpl::is_void_<CONCRETE>,
                     AbstractAssembler<ELEMENT_DIM, SPACE_DIM, 1u>,
                     typename AssemblerTraits<CONCRETE>::CMT_CLASS>::type
            INTERPOLATE_CLASS;
};

New SimpleLinearEllipticSolver 


#include "AbstractLinearEllipticPde.hpp"
// the 2 new parents new to be included 
#include "AbstractAssemblerSolverHybrid.hpp"
#include "AbstractStaticLinearPdeSolver.hpp"  // or AbstractDynamicinearPdeSolver


template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
class SimpleLinearEllipticSolver 
     // new parents: See assembler/solver description page for what NORMAL means.
     // Second should be either AbstractStaticLinearPdeSolver or AbstractDynamicLinearPdeSolver
   : public AbstractAssemblerSolverHybrid<ELEMENT_DIM, SPACE_DIM, 1, NORMAL>,  
     public AbstractStaticLinearPdeSolver<ELEMENT_DIM, SPACE_DIM, 1>     
{
protected:
    AbstractLinearEllipticPde<ELEMENT_DIM,SPACE_DIM>* mpEllipticPde;

    
    virtual c_matrix<double, 1*(ELEMENT_DIM+1), 1*(ELEMENT_DIM+1)> ComputeMatrixTerm(
        c_vector<double, ELEMENT_DIM+1>& rPhi,
        c_matrix<double, SPACE_DIM, ELEMENT_DIM+1>& rGradPhi,
        ChastePoint<SPACE_DIM>& rX,
        c_vector<double,1>& rU,
        c_matrix<double,1,SPACE_DIM>& rGradU,
        Element<ELEMENT_DIM,SPACE_DIM>* pElement)
    {
       .. 
    }

    
    virtual c_vector<double,1*(ELEMENT_DIM+1)> ComputeVectorTerm(
        c_vector<double, ELEMENT_DIM+1>& rPhi,
        c_matrix<double, SPACE_DIM, ELEMENT_DIM+1>& rGradPhi,
        ChastePoint<SPACE_DIM>& rX,
        c_vector<double,1>& rU,
        c_matrix<double,1,SPACE_DIM>& rGradU,
        Element<ELEMENT_DIM,SPACE_DIM>* pElement)
    {
       ..
    }
        

    virtual c_vector<double, ELEMENT_DIM> ComputeVectorSurfaceTerm(const BoundaryElement<ELEMENT_DIM-1,SPACE_DIM>& rSurfaceElement,
            c_vector<double, ELEMENT_DIM>& rPhi,
            ChastePoint<SPACE_DIM>& rX)
    {
       ..
    }

    //  **** ADD THE FOLLOWING METHOD ****. 
    // It implements a pure method on the solver parent class by
    // calling a method on the hyrid parent.
    
    /**
     * Delegate to AbstractAssemblerSolverHybrid::SetupGivenLinearSystem.
     *  @param currentSolution The current solution which can be used in setting up
     *   the linear system if needed (NULL if there isn't a current solution)
     *  @param computeMatrix Whether to compute the LHS matrix of the linear system
     *   (mainly for dynamic solves).
     */
    void SetupLinearSystem(Vec currentSolution, bool computeMatrix)
    {
        SetupGivenLinearSystem(currentSolution, computeMatrix, this->mpLinearSystem);
    }
    
public:

    // Constuctor calls base class constructors
    SimpleLinearEllipticSolver(AbstractTetrahedralMesh<ELEMENT_DIM,SPACE_DIM>* pMesh,
                               AbstractLinearEllipticPde<ELEMENT_DIM,SPACE_DIM>* pPde,
                               BoundaryConditionsContainer<ELEMENT_DIM,SPACE_DIM,1>* pBoundaryConditions,
                               unsigned numQuadPoints = 2)
        : AbstractAssemblerSolverHybrid<ELEMENT_DIM,SPACE_DIM,1,NORMAL>(pMesh,pBoundaryConditions,numQuadPoints),
          AbstractStaticLinearPdeSolver<ELEMENT_DIM,SPACE_DIM,1>(pMesh)
    {
       mpEllipticPde = pPde;
    }


    void InitialiseForSolve(Vec initialSolution = NULL)
    {
       // Only change is the 'AbstractLinearPdeSolver<ELEMENT_DIM,SPACE_DIM,1>::'
       AbstractLinearPdeSolver<ELEMENT_DIM,SPACE_DIM,1>::InitialiseForSolve(initialSolution);
       
       assert(this->mpLinearSystem);
       this->mpLinearSystem->SetMatrixIsSymmetric(true);
       this->mpLinearSystem->SetKspType("cg");
    }
};

Attachments