SystemNSE3D.h

#ifndef __SYSTEMNSE3D__
#define __SYSTEMNSE3D__

#include <SquareMatrix3D.h>
#include <FEVectFunct3D.h>
#include <NSE_MultiGrid.h>
#include <ItMethod.h>
#include <AssembleMat3D.h>

#ifdef _MPI
//#include "mpi.h"
#include <ParFEMapper3D.h>
#include <ParFECommunicator3D.h>
#include <ParDirectSolver.h>
#endif

#ifdef _OMPONLY
#include <ParDirectSolver.h>
#endif

class TSystemNSE3D
{
  protected:

#ifdef _MPI
    TParFEMapper3D **ParMapper_U, **ParMapper_P;
    TParFECommunicator3D **ParComm_U, **ParComm_P;
    MPI_Comm Comm;
    TParDirectSolver *DS;
#endif
    
#ifdef _OMPONLY
    TParDirectSolver *DS;
#endif
    
    int N_Levels;
    
    int Start_Level;  
    
    int N_TotalDOF, N_U, N_P, N_Active, N_DirichletDof;
    
    int *N_FreeSurfFaces;
    
    int **FreeSurfCellNumbers, **FreeSurfJointNumbers;
    
    TFESpace3D **U_Space, **P_Space, *FeSpaces[5];
    
    TFEVectFunct3D **Velocity;
    
    TFEFunction3D **Pressure;    

    double **SolArray, **RhsArray, *RHSs[4];

    TFESpace3D *fesp[3], *fesprhs[3], *fesp_aux[3];
    TFEFunction3D *fefct[7], **fefct_aux;
  
    // Assembling */
    TAssembleMat3D **AMatRhsAssemble, **AMatAssembleNonLinear;
    
    int Disctype;
       
    int NSEType; 
           
    CoeffFct3D *LinCoeffs[1];    
    
    TAuxParam3D **NSEaux, **NSEaux_error;
    
    MatVecProc *MatVect;  
    
    DefectProc *Defect;   
    
    int SOLVER;
       
    int N_Matrices;

    TSquareStructure3D **sqstructureA;

    TStructure3D **structureB, **structureBT;
    
    TSquareMatrix3D **SqmatrixA11, **SqmatrixA12, **SqmatrixA13;
    TSquareMatrix3D **SqmatrixA21, **SqmatrixA22, **SqmatrixA23;
    TSquareMatrix3D **SqmatrixA31, **SqmatrixA32, **SqmatrixA33, *SQMATRICES[15];
    TSquareMatrix **sqmatrices;
    TSquareMatrix3D **SqmatrixF11, **SqmatrixF22, **SqmatrixF33;
    
    TMatrix3D **MatrixB1, **MatrixB2, **MatrixB3, **MatrixB1T, **MatrixB2T, **MatrixB3T, *MATRICES[12];
    TMatrix **matrices;
    
    BoundCondFunct3D *BoundaryConditions[3];
  
    BoundValueFunct3D *BoundaryValues[3];
        
    TDiscreteForm3D *DiscreteFormARhs, *DiscreteFormNL;

    int N_aux;
    double Parameters[2], *Itmethod_sol, *Itmethod_rhs;
    TNSE_MultiGrid *MG;
    TNSE_MGLevel *MGLevel;
    TItMethod *Itmethod, *prec;
 
    
  private:
   void UpdateUpwind(int i);
    
//    void UpdateLPS();
   
  public:
     TSystemNSE3D(int N_levels, TFESpace3D **velocity_fespace, TFESpace3D **presssure_fespace, TFEVectFunct3D **velocity, 
                     TFEFunction3D **pressure, double **sol, double **rhs, int disctype, int nsetype, int solver);

//     ~TSystemMatNSE3D();

    void Init(CoeffFct3D *lincoeffs, BoundCondFunct3D *BoundCond, BoundValueFunct3D *U1BoundValue, 
              BoundValueFunct3D *U2BoundValue, BoundValueFunct3D *U3BoundValue);
    
    void Assemble();
    
    void AssembleNonLinear(double **sol, double **rhs);
    
    void GetResidual(double *sol, double *rhs, double *res, double &impuls_residual, double &residual);
    
    void  Solve(double *sol, double *rhs);
  
    void MeasureErrors(DoubleFunct3D *ExactU1, DoubleFunct3D *ExactU2, DoubleFunct3D *ExactU3, DoubleFunct3D *ExactP,
                       double *u_error, double *p_error);
    
    void FindFreeSurfJoints(int level, int Region_ID);
};

#endif