MGsystem.h

//===========================================================================
// SINTEF LSMG library - version 1.1
//
// Copyright (C) 2000-2005 SINTEF ICT, Applied Mathematics, Norway.
//
// This program is free software; you can redistribute it and/or          
// modify it under the terms of the GNU General Public License            
// as published by the Free Software Foundation version 2 of the License. 
//
// This program is distributed in the hope that it will be useful,        
// but WITHOUT ANY WARRANTY; without even the implied warranty of         
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the          
// GNU General Public License for more details.                           
//
// You should have received a copy of the GNU General Public License      
// along with this program; if not, write to the Free Software            
// Foundation, Inc.,                                                      
// 59 Temple Place - Suite 330,                                           
// Boston, MA  02111-1307, USA.                                           
//
// Contact information: e-mail: tor.dokken@sintef.no                      
// SINTEF ICT, Department of Applied Mathematics,                         
// P.O. Box 124 Blindern,                                                 
// 0314 Oslo, Norway.                                                     
//
// Other licenses are also available for this software, notably licenses
// for:
// - Building commercial software.                                        
// - Building software whose source code you wish to keep private.        
//===========================================================================
#ifndef _MGSYSTEM_H_
#define _MGSYSTEM_H_

#include <LSsystem.h>
#include <UCBsplineSurface.h>

class MGsystem {

    // Common temporaries for both FMG and V-cycle (but not for solveAscend)
    std::vector<double> domain_;    // defailt is umin, vmin, umax, vmax (see also LSsystem)
    std::vector<LSsystem*> systems_;
    GenMatrix<UCBspl_real> r_;
    GenMatrix<UCBspl_real> e_refined_;
  
    // Array of solution vectors (pointers)
    typedef boost::shared_ptr<GenMatrix<UCBspl_real> > shared_gm;
    std::vector<shared_gm> e_arr_;

    int noIter_;
    int noIterCoarsest_;

    double lambdaFac_;

    int m0_, n0_; // given to init multilevel 
    void calculateStartVector() const;

    double umin() const {return domain_[0];}
    double vmin() const {return domain_[1];}
    double umax() const {return domain_[2];}
    double vmax() const {return domain_[3];}

    void reserveResidual(); // assumes that init has been done

    void residual(const LSsystem& lssys, GenMatrix<UCBspl_real>& r) const;
    void residual_restrictionInjection(const LSsystem& lssys, GenMatrix<UCBspl_real>& r) const;
    void residual_restrictionFullWeightingLaplace(const LSsystem& lssys, 
                                                  GenMatrix<UCBspl_real>& r) const;
    void residual_restrictionSplines(const LSsystem& lssys, GenMatrix<UCBspl_real>& r) const;

    // Can also be used with a start vector (e.g. to adjust smoothness)
    // ???? Is not using private members now
    void relaxAndCorrectVcycle(int noLevelsInVcycle);

    // Used to control allocation of right hand side (RHS) at the finest level.
    // This is used to avoid building the RHS again if a solveVcycle or relax is
    // run afterwards to improve the solution (or adjust smoothness)
    // false = initial status before any solvers have been run
    // true  = the RHS at the finest level has been calculated
    bool rhsFinestLevel_;

public:
    MGsystem() {lambdaFac_=1.0; noIter_=10; noIterCoarsest_=50; rhsFinestLevel_=false;}
    ~MGsystem() {cleanup();}
  
    typedef boost::shared_ptr<std::map<int, double> > shared_map;
    void initMultilevel(int m0, int n0, int h,
                        boost::shared_ptr<std::vector<double> > U,
                        boost::shared_ptr<std::vector<double> > V,
                        boost::shared_ptr<std::vector<double> > Z,
                        boost::shared_ptr<GenMatrix<UCBspl_real> > x, // solution vector
                        shared_map weights = shared_map(new std::map<int, double>) );
    //boost::shared_ptr<std::map<int, double> > weights = boost::shared_ptr<std::map<int, double> >(new std::map<int, double>) );

  
              
    void solveFMG();
  
    void solveAscend();

    //void solveAscend(double err);

    void solveVcycle();

    void relax(int numberOfIterations, int relaxType = 3);

    UCBspl::SplineSurface getSplineSurface() const;

    double residual_l2(bool scaled=false) const;
    double residual_linf() const;

    void setDomain(double umin, double vmin, double umax, double vmax);

    void getDomain(double& umin, double& vmin, double& umax, double& vmax);

    void setSmoothingFactor(double smoothingFac);

    bool addPoint(double u, double v, double z, double weight=1.0);

    void setWeight(int pointIndex, double weight);
  
    void setNumberOfIterations(int noIter) {noIter_=noIter;}

    void setNumberOfIterationsCoarsest(int noIter) {noIterCoarsest_ = noIter;}
  
    void cleanup();

    /* ??? Initialize solution vector with \e value */
    void qweInitX(double value) {
        systems_[systems_.size()-1]->x_->fill(value);
    }

    // ??? temporary function for use in LSsystem::CGPrecond
    /* Set right hand side in equation system. */
    void qweSetRHS(const GenMatrix<UCBspl_real>& t_k) {
        systems_[systems_.size()-1]->b_.init(t_k);
        rhsFinestLevel_ = true;
    }

};

#endif

---