SystemNSE2D.h

#ifndef __SYSTEMNSE2D__
#define __SYSTEMNSE2D__

#include <SquareMatrix2D.h>
#include <AssembleMat2D.h>

class TSystemNSE2D
{
  protected:

    int N_U, N_P, N_Active, N_DirichletDof;
    
   double *Sol, *Rhs, *RHSs[3];   
    
    
    TFESpace2D *FeSpaces[5];
    
    TFEFunction2D *FeFct[5];    
    
    int Disctype;
       
    int NSEType; 
           
    CoeffFct2D *LinCoeffs[1];    
    
    TAuxParam2D *NSEaux, *NSEaux_error;
    
    DefectProc *Defect;   
    
    int Solver;
       
    int N_Matrices;

    TSquareStructure2D *sqstructureA;

    TStructure2D *structureB, *structureBT;
    
    TSquareMatrix2D *SqmatrixA11, *SqmatrixA12, *SqmatrixA21, *SqmatrixA22, *SQMATRICES[9];
    TSquareMatrix **sqmatrices;
  
    TMatrix2D *MatrixB1, *MatrixB2, *MatrixB1T, *MatrixB2T, *MATRICES[8];
    TMatrix **matrices;
    
    BoundCondFunct2D *BoundaryConditions[2];
  
    BoundValueFunct2D *BoundaryValues[2];
        
    TDiscreteForm2D *DiscreteFormARhs, *DiscreteFormNL, *DiscreteFormRhs;

    TAssembleMat2D *AMatRhsAssemble, *AMatAssembleNonLinear;  
   
   
#ifdef __PRIVATE__  

    TSquareStructure2D *sqstructureL;

    TStructure2D *structure_G, *structure_tilde_G;      
    
    TSquareMatrix2D *MatricesL;
    
    TMatrix2D *Matrices_tilde_G11, *Matrices_tilde_G22, *Matrices_G11, *Matrices_G22;  
#endif 
 
  private:
   void UpdateUpwind();
    
   void UpdateLPS();
  public:
     TSystemNSE2D(TFESpace2D *velocity_fespace, TFESpace2D *presssure_fespace, TFEVectFunct2D *Velocity, 
                  TFEFunction2D *p, double *sol, double *rhs, int disctype, int nsetype, int solver
#ifdef __PRIVATE__  
                  ,TFESpace2D *Projection_space, TFESpace2D *Stress_FeSpace, TFESpace2D *Deformation_FeSpace
#endif         
                );

    ~TSystemNSE2D();

    void Init(CoeffFct2D *lincoeffs, BoundCondFunct2D *BoundCond, BoundValueFunct2D *U1BoundValue, BoundValueFunct2D *U2BoundValue,
              TAuxParam2D *aux, TAuxParam2D *auxerror) ;
    
    
    void Assemble(double *sol, double *rhs);
    
    void AssembleNonLinear(double *sol, double *rhs);
    
    void AssembleRhsOnly(double *sol, double *rhs);
    
    void GetResidual(double *sol, double *rhs, double *res);
    
    void  Solve(double *sol, double *rhs);
  
    void MeasureErrors(DoubleFunct2D *ExactU1, DoubleFunct2D *ExactU2, DoubleFunct2D *ExactP,
                                    double *u_error, double *p_error);
};

#endif