SystemTCD3D.h

#ifndef __SYSTEMTCD3D__
#define __SYSTEMTCD3D__

#include <SquareMatrix3D.h>
#include <SystemCD3D.h>

class TSystemTCD3D : public TSystemCD3D
{
  protected:
#ifdef _MPI
    TParDirectSolver *TDS;
#endif   

    TSquareMatrix3D **sqmatrixM;
    
    double *B;   
   
    double *defect;   
    
    double gamma;   
    
    bool factorize;

    
    TAssembleMat3D **MMatRhsAssemble;
    
//     /** Stiffness part of the SUPG matrix */
//     TSquareMatrix3D *sqmatrixK;    
//     
//     /** time-consistent part of the SUPG matrix */
//     TSquareMatrix3D *sqmatrixS;
    
    TDiscreteForm3D *DiscreteFormMRhs, *DiscreteFormRhs; 
    
    bool SystMatAssembled;
    
  public:
     TSystemTCD3D(int N_levels, TFESpace3D **fespaces, double **sol, double **rhs, int disctype, int solver);

    ~TSystemTCD3D();

    void Init(CoeffFct3D *BilinearCoeffs, BoundCondFunct3D *BoundCond, BoundValueFunct3D *BoundValue, TAuxParam3D *aux);
    
    TSquareMatrix3D **GetAMatrix()
    { return sqmatrixA; }
    
    void AssembleMRhs(); 
    
    void AssembleARhs();   
    
//     /** M = M + (tau*THETA1)*A */ 
//     /** B = (tau*THETA1)*rhs +(tau*THETA2)*oldrhs + [ M - (tau*THETA2)A]*oldsol */  
    void AssembleSystMat(double *oldrhs, double *oldsol, double *rhs, double *sol
#ifdef _MPI
                         , double **Rhs_array
#endif
                         );
    
    void RestoreMassMat();
    
    void Solve(double *sol);  
    
    double GetResidual(double *sol);
    
    double value(double *sol,int N){
      int i;
      double sum=0.;
      for(i=0;i<N;i++)  sum+=sol[i];
      return sum;
    }
    
};

#endif