#-*- coding: utf-8 -*-
__doc__ = """

NAVIERSTOKES
=============

The main class for the navier-stokes method. 

* Note: No Turbulence scheme implemented (only laminar flow)

Description
-----------
This class wraps the naviers-stokes solver for time-stepping with a given
dirichlet velocity boundary condition. However, the solver is also capable of 
solving the problem with a pressure-outlet (neumann b.c.) if the dirichlet 
velocity boundary condition at the exit is unknown.

The wrapped NS solver is Finite Element Method Solver for incompressible 
**laminar** fluid. This means that no turbulence scheme is implemented. The

Methodology
-----------
1. Define the NS problem by inputing the parameters of the geometry, the fluid,
   the solver and the probe mesh. 

2. Retrive the boundary coordinates of the FE mesh.

3. Determine the dirichlet boundary condition at those node.

4. Use the dirichlet b.c, new mesh position (and velocity) to evolve the problem.

5. [hybrid] Retrive the vorticity data from the probes.   

6. Repeat.  

Implemented NS Algorithms
-------------------------
- Chorin projection scheme (i.e Non-Incremental pressure correction)
- Incremental pressure correction scheme (ipcs)

:First Added:   2013-12-10
:Last Modified: 2014-03-06                       
:Copyright:     Copyright (C) 2013 Lento Manickathan **pHyFlow**
:License:       GNU GPL version 3 or any later version
"""

#   Copyright (C) 2013 Lento Manickathan                                                                         
#   
#   This file is part of pHyFlow.                                                                                                      
#   
#   pHyFlow is free software: you can redistribute it and/or modify                                                                    
#   it under the terms of the GNU Lesser General Public License as published by                                                       
#   the Free Software Foundation, either version 3 of the License, or                                                                 
#   (at your option) any later version.                                                                                               
#   
#   pHyFlow 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 Lesser General Public License for more details.                                                                               
#      
#   You should have received a copy of the GNU Lesser General Public License                                                          
#   along with pHyFlow. If not, see <http://www.gnu.org/licenses/>.    

__all__ = ['EulerianSolver']

# External packages
import numpy as _numpy
import fenicstools as _fenicstools

# Import pHyFlow packages
from pHyFlow.aux.customDecorators import simpleGetProperty

# import eulerian packages
from pHyFlow.eulerian import base as _base
from pHyFlow.eulerian import eulerianOptions

class EulerianSolver(object):
    r"""
    Class containing all the overhead functions and managing the eulerian
    solver using the FEniCS/DOLFIN Finite Element Solver.
    
    Usage
    -----
    .. code-block:: python
    
        eulerian = EulerianSolver(geometry = {'mesh': meshFile,
                                              'boundaryDomains': boundaryDomainsFile,
                                              'cmGlobal': cmGlobal,
                                              'thetaLocal': thetaLocal},
                                  probeGrid = {'origin': origin,
                                               'L': L
                                               'N': N,
                                               'circulationProbes': xyCirculationProbeCoordinates},
                                  uMax=uMax,nu=nu,cfl=cfl,deltaT,
                                  solverParams={'solver':'ipcs')
        
    Parameters
    ----------
    geometry : dict    
               dictionary containing {'mesh','boundaryDomains','cmGlobal',
               'thetaLocal'} parameters of the geometry.
               
               'mesh' : str
                        the mesh data filename (.xml.gz) generated by a
                        FE mesh generateor (GMSH) and converted to XML format.
                      
               'boundaryDomains' : str
                                   the facet boundary domain mesh data file 
                                   location. The facet boundary file should be
                                   marked (int format) according to 
                                   .. py:module:: pHyFlow.navierStokes.nsOptions
                                   as shown below:
                                 
                                   (1) : Fluid Domain
                                   (2) : No-Slip boundary
                                   (3) : External dirichlet boundary
                                   (4) : [optional] Pressure outlet.

                                   * Note: Pressure outlet is optional.
                                 
               'cmGlobal' : numpy.ndarray(float64), shape (2,)
                            the :math:`x,y` position of the mesh local reference
                            point (0.,0.) in the global coordinates.
                          
               'thetaLocal' : float, unit (rad)
                              the local rotational angle :math:`\theta` of the 
                              mesh domain. Therefore, the rotation will be done
                              about local reference point (0.,0.), i.e cmGlobal 
                              in the global coordinate system.
                            
    probeGrid : dict
                dictionary containing all the parameters of the probe grid for
                extracting the vorticity data.                             

                'origin' : numpy.ndarray(float64), shape (2,)
                           the :math:`x,y` coordinate of the origin of the probe 
                           mesh in the **local** coordinate system of the fluid
                           mesh.
                         
                'L' : numpy.ndarray(float64), shape (2,)
                      the width :math:`L_x` and the height :math:`L_y` of the 
                      probe mesh.
                
                'N' : numpy.ndarray(float64), shape (2,)
                      the number of probes :math:`N_x,N_y` in the :math:`x,y` 
                      direction.
                      
                'circulationProbes': numpy.ndarray(float64), shape(2,nProbes), orientation (anti-clockwise)
                                     the :math:`x,y` local coordinates of the
                                     polygon defining the area where the 
                                     circulation is to be calculated. The
                                     circulation of the area defined by the 
                                     polygon is calculated using stokes theorem
                                     of circulation. Therefore, should be 
                                     in **ANTI-CLOCKWISE**.
                                     
                                     * Note: for hybrid coupling, the circulation
                                     probes should be located at the same location
                                     as the outter boundary of the interpolation
                                     domain.
    
    uMax : float
           the maximum fluid velocity :math:`U_{max}`.
        
    nu : float
         the fluid kinematic viscosity :math:`\nu`.
        
    cfl : float
          the :math:`CFL` stability parameter. For explicit time-stepping,
          :math:`CFL\le1`
        
    solverParams : dict
                   dictionary file containing all the solver parameters.
                   
                   'solver' : str
                              type of NS solver that should used to solve
                              the incompressible laminar problem. Available
                              solver can be found in .. py:module:: pHyFlow.navierStokes.nsOptions:
                              
                              'ipcs' : [default] Incremental pressure-correction scheme
                              'chorin': Chorin projection scheme (non-incremental)
        
    
    Attribute
    ---------
    deltaT
    deltaTMax
    cfl
    cmGlobal
    hMin    
    nu
    probeGridMesh
    probeGridParams
    solverParams
    t
    thetaLocal
    tStep
    uMax
        
    __boundaryDomainsFile : str
                            the facet boundary domain mesh data file location.                 
    
    __cfl : float
            the :math:`CFL` stability parameter.

    __cmGlobal : numpy.ndarray(float64), shape (2,)
                 the :math:`x,y` position of the mesh local reference point 
                 (0.,0.) in the global coordinates.
                 
    __meshFile : str
                 the mesh data filename (.xml.gz) generated by a FE mesh 
                 generateor (GMSH) and converted to XML format.
                        
    __nu : float
           the fluid kinematic viscosity :math:`\nu`.
    
    __probeGrid : dict
                dictionary containing all the parameters of the probe grid for
                extracting the vorticity data.                             

                'origin' : numpy.ndarray(float64), shape (2,)
                           the :math:`x,y` coordinate of the origin of the probe 
                           mesh in the **local** coordinate system of the fluid
                           mesh.
                         
                'L' : numpy.ndarray(float64), shape (2,)
                      the width :math:`L_x` and the height :math:`L_y` of the 
                      probe mesh.
                
                'N' : numpy.ndarray(float64), shape (2,)
                      the number of probes :math:`N_x,N_y` in the :math:`x,y` 
                      direction.
                         
    __probes : fenicstools.Probes.Probes
               the `probe` class containing all the information to `eval`
               the vorticity function at the probe grid coordinates.
    
    __solver : pHyFlow.navierStokes.base.solverBase
               the NS solver class containing all the information to solver
               the incompressible laminar navier-stoke problems to the given
               input geometry, fluid and solver parameters.

    __solverParams : dict
                     dictionary file containing all the solver parameters.
                     'solver' : str
                                type of NS solver that should used to solve
                                the incompressible laminar problem.
    
    __thetaLocal : float
                   the local rotational angle :math:`\theta` of the mesh domain.
                   Therefore, the rotation will be done about local reference 
                   point (0.,0.), i.e cmGlobal in the global coordinate system.

    __uMax : float
             the maximum fluid velocity :math:`U_{max}`.

    __xyProbes : numpy.ndarray(float64), shape (__probeGrid_N[0], __probeGrid_N[1])
                 the mesh grid of probes generted by the parameters of the
                 probe. The probe mesh grid origin at **__probeGrid_origin**,
                 with length **__probeGrid_L** with **__probeGrid_N** probes
                 is :math:`x,y` direction.
    
    Methods
    -------
    getCoordinates
    setVelocity

    getBoundaryCoordinates
    
    evolve

    getVorticity

    getProbeGrid
    
   

    :First Added:   2013-12-10
    :Last Modified: 2014-02-21                         
    :Copyright:     Copyright (C) 2013 Lento Manickathan **pHyFlow**
    :License:       GNU GPL version 3 or any later version   
   
    """
    """
    Revisions
    ---------
    2013-12-10, Lento Manickathan
        - First added
        
    2014-02-21, Lento Manickathan
        - Documentation update
        - new boundary search algorithm
        - thetaGlobal changed to thetaLocal

    """
    
    def __init__(self,geometry, probeGrid, uMax, nu, cfl,deltaT='max',
                 solverParams={'solver':eulerianOptions.SOLVER['default']}):
        
        #---------------------------------------------------------------------
        # Check/Set input parameters
            
        # Initialize the parameters
        self.__set_variables('geometry', geometry)
        self.__set_variables('probeGrid', probeGrid)
        self.__set_variables('uMax', uMax)
        self.__set_variables('nu',nu)            
        self.__set_variables('cfl',cfl)  
        self.__set_variables('solverParams', solverParams)
        #---------------------------------------------------------------------
        
        #---------------------------------------------------------------------
        # Choose the solver: (Chorin or other solvers)
        # - Various Navier-stokes solving algorithms can be imported. All the
        # algorithms should be stored in './solvers'. The solver is then
        # imported in as part of this module into -> 'self.solver'.
        if self.__solverParams['solver'] == 'chorin':
            self.__solver = _base.chorin(self.__meshFile,self.__boundaryDomainsFile,
                                       self.__nu,self.__cfl,self.__uMax)
        elif self.__solverParams['solver'] == 'ipcs':
            self.__solver = _base.ipcs(self.__meshFile,self.__boundaryDomainsFile,
                                       self.__nu,self.__cfl,self.__uMax)
        else:
            raise NotImplementedError('Solver type not implemented or unknown !!')            
        #---------------------------------------------------------------------

        #---------------------------------------------------------------------            
        # Define the probe grid, w.r.t to the local mesh
        # * Note: the probes are defined/generated in the local coordinates system.
        #         Therefore, we must probe the mesh, before we update the mesh
        #         to global coordinate system.
        self.__xyProbes = _numpy.meshgrid(_numpy.linspace(0,self.__probeGrid['L'][0],
                                                          self.__probeGrid['N'][0]) + self.__probeGrid['origin'][0],
                                          _numpy.linspace(0,self.__probeGrid['L'][1],
                                                          self.__probeGrid['N'][1]) + self.__probeGrid['origin'][1])
                     
        # Initialize the probes                     
        self.__probes = _fenicstools.Probes(_numpy.append(self.__xyProbes[0].reshape(-1,1),
                                                          self.__xyProbes[1].reshape(-1,1),axis=1).flatten(),
                                            self.__solver.X)
                                            
        # Initialize the circulation probes
        # Determine the center points.
        circulationProbeEdgeCoordinates = self.__probeGrid['circulationProbes']                                            
        circulationProbesMidCoordinates = 0.5*(circulationProbeEdgeCoordinates[:,1:]+circulationProbeEdgeCoordinates[:,:-1])
                                            
  
        # Determine the circulation probes parameters
        self.__circulationProbes_dl = circulationProbeEdgeCoordinates[:,1:] \
                                        - circulationProbeEdgeCoordinates[:,:-1]
                         
        # Initialize the circulation probes                             
        self.__circulationProbes = _fenicstools.Probes(circulationProbesMidCoordinates.T.flatten(),
                                                       self.__solver.V)
                                                              
        #---------------------------------------------------------------------
                                            
        #---------------------------------------------------------------------                                            
                                            
        #---------------------------------------------------------------------
        # Update mesh position (new rotational angle and global position)

        # The rotating is done around the local reference point (0.,0.)
        #dThetaLocal = self.__thetaLocal - 0.0 # new - old angle
        #self.__solver.rotateMesh(dThetaLocal) # rotate by dThetaLocal
        self.__solver.updateMeshPosition(self.__cmGlobal, self.__thetaLocal)
        #---------------------------------------------------------------------

        #---------------------------------------------------------------------                                            
        # Define and initialize the time step counters and absolute time

        # assign the delta T        
        self.__set_variables('deltaT',deltaT)

        # define the time step counter
        self.__tStep = 0
        
        # define the current time
        self.__t = 0.0
                          
        # We need to calculate the vorticity  
        self.__solver.recalculateVorticityFlag = True
        self.__reprobeVorticityFlag = True 
        
        # plotting flags
        self.plotVelocity = False
        self.plotPressure = False
        self.plotVorticity = True # most informative data
        #---------------------------------------------------------------------
                                    
                                    
    def evolve(self,vxBoundary,vyBoundary,cmGlobalNew,thetaLocalNew,cmDotGlobal,
               thetaDotLocal):
        r"""
        Function to evolve the Navier-Stokes by one step with the :math:`x,y`
        velocity boundary conditions at the navier-stokes dirichlet boundary.
        
        The dirichlet boundary is identified by 
        .. py:module:: pHyFlow.navierStokes.nsOptions.ID_EXTERNAL_BOUNDARY
        in the **boundaryDomains**. The coordinate of the dirichlet boundary
        can be obtained from .. py:module:: getBoundaryCoordinates
        
        The function will calculate the new velocity and the pressure fields.
        
        * Note : [not implemented] Moving mesh.
            
        
        Usage
        -----
        .. code-block:: python
        
            evolve(vxBoundary,vyBoundary,cmGlobal,thetaLocal,cmDotGlobal,
                   thetaDotLocal)
        
        Parameters
        ----------
        vxBoundary : numpy.ndarray(float64), shape (nDOFs,)
                     the :math:`x` component of the dirichlet velocity boundary
                     condition at the navierstokes DOF boundaries.

        vyBoundary : numpy.ndarray(float64), shape (nDOFs,)
                     the :math:`y` component of the dirichlet velocity boundary
                     condition at the navierstokes DOF boundaries.                                  
        
        cmGlobalNew : numpy.ndarray(float64), shape (2,)
                      the new :math:`x,y` mesh position of the navierStokes domain.
                      * Note: Note implemented yet. Does not update the position.
        
        thetaLocalNew : float
                        the new mesh rotational angle :math:`\theta`.
                        * Note: Note implemented yet. Does not update the position.
        
        cmDotGlobal : numpy.ndarray(float64), shape (2,)
                      the current :math:`x,y` mesh velocity.
                      * Note: Note implemented yet. Does not update the position.
        
        thetaDotGlobal : float
                         the current mesh rotational velocity :math:`\dot{\theta}`. 
                         * Note: Note implemented yet. Does not update the position.
                         
        Returns
        -------
        None.
        
        Attribute
        ---------
        __solver : 
                    time steps the solver and updates the parameters of the
                    solver.
                    
                    __solver.bcVelocity : list
                                          update the velocity b.c
                                            
                    __solver.u0, __solver.u1 : dolfin.function
                                               update the new and old velocity
                                               function.
                    __solver.p0, __solver.p1 : dolfin.function
                                               update the new and old pressure
                                               function.
                
        
        :First Added:   2013-12-19
        :Last Modified: 2014-02-21
        :Copyright:     Copyright (C) 2014 Lento Manickathan **pHyFlow**
        :Licence:       GNU GPL version 3 or any later version 
        
        """
        # update the mesh position [moving boundary problem not implemented !!]
        #self.__thetaLocal = thetaLocalNew
        #self.__cmGlobal = cmGlobalNew
        #self.__solver.updateMeshPosition(cmGlobalNew, thetaLocalNew)
        
        # Set boundary conditions
        self.__solver.boundaryConditions(vxBoundary,vyBoundary)
        
        # Solve the navier-stokes problem
        self.__solver.solve()
        
        # update the time
        self.__advanceTime()  
   
   
    def Forces(self):
        r"""
        Function to calculate the Pressure gradient forces (N).
        
        If the free-stream flow is in x-direction (1,0):
            
            Fx : Total drag force
            Fy : Total lift force
        
        Usage
        -----
        .. code-block:: python
        
            F = Forces()
            
        Parameters
        ----------
        None
        
        Returns
        -------
        F : numpy.ndarray(float64), shape (2,), unit (N)
            the :math:`x,y` component of the total forces acting on the
            geometry (defined by the no-slip boundary).
            
        :First Added:   2014-03-05
        :Last Modified: 2014-03-05
        :Copyright:     Copyright (C) 2014 Lento Manickathan **pHyFlow**
        :Licence:       GNU GPL version 3 or any later version  
        """
        return self.__solver.Forces()  
        

    def FrictionalForces(self):
        r"""
        Function to calculate the frictional forces. 
        
        If the free-stream flow is in x-direction (1,0):
            
            Ffx : Friction drag
            Ffy : Friction induced lift
        
        Usage
        -----
        .. code-block:: python
        
            Ff = FrictionalForces()
            
        Parameters
        ----------
        None
        
        Returns
        -------
        Ff : numpy.ndarray(float64), shape (2,), unit (N)
             the :math:`x,y` component of the frictional forces acting on the
             geometry (defined by the no-slip boundary).
            
        :First Added:   2014-03-05
        :Last Modified: 2014-03-05
        :Copyright:     Copyright (C) 2014 Lento Manickathan **pHyFlow**
        :Licence:       GNU GPL version 3 or any later version         
        
        """
        return self.__solver.FrictionalForces()
  

    def getBoundaryCoordinates(self):
        """
        Returns the *global* boundary DOF coordinates :math:`x,y_{boundary}` of 
        the vector function space :math:`\mathbf{V}`, where the dirichlet
        velocity boundary condition should be applied.
        
        The boundary coordinates is located with the help of the 
        **boundaryDomains** mesh function. The dirichlet velocity boundary
        DOF coordinates is found by a internal FEniCS search initialization
        of the module.
        
        The Dirichlet boundary is ID is found in 
        .. py:module:: pHyFlow.navierStokes.nsOptions.ID_EXTERNAL_BOUNDARY
        
            External boundary: 3
        
        Usage
        -----
        .. code-block:: python
        
            xBoundary, yBoundary = getBoundaryCoordinates()
            
        or
        
        .. code-block :: python
        
            xyBoundary = getBoundaryCoordinates()
            
        Parameters
        ----------
        None
            
        Returns
        -------
        xyBoundary : numpy.ndarray(float64), shape (2,nDOFs)
                     the global :math:`x,y` coordinates of the dirichlet boundary
                      where the external velocity is applied.
        
        Attributes
        ----------
        None changed.
        
        :First Added:   2013-12-18
        :Last Modified: 2014-03-06
        :Copyright:     Copyright (C) 2013 Lento Manickathan **pHyFlow**
        :Licence:       GNU GPL version 3 or any later version    
        """
        
        # return the global coordinates of the vector dof boundaries
        return self.__solver.vectorDOF_boundaryCoordinates()
        
 
    def getBoundaryVelocity(self):
        """
        Returns the boundary velocity of the dirichlet velocity boundary
        domain.
        
        The Dirichlet boundary is ID is found in 
        .. py:module:: pHyFlow.navierStokes.nsOptions.ID_EXTERNAL_BOUNDARY
        
            External boundary: 3
        
        Usage
        -----
        .. code-block:: python
        
            vxBoundary, vyBoundary = getBoundaryVelocitys()
            
        Parameters
        ----------
        None
            
        Returns
        -------
        None returned.
        
        Attributes
        ----------
        None changed.
        
        :First Added:   2013-02-28
        :Last Modified: 2014-02-28
        :Copyright:     Copyright (C) 2014 Lento Manickathan **pHyFlow**
        :Licence:       GNU GPL version 3 or any later version    
        """
        
        # Get the boundary coordinates of the mesh (local coordinate system)
        boundary_VectorDOFIndex = self.__solver.boundary_VectorDOFIndex
        
        # in the global coordinate system
        vx = self.__solver.u1.vector()[boundary_VectorDOFIndex[0]].copy()
        vy = self.__solver.u1.vector()[boundary_VectorDOFIndex[1]].copy()
        
        # return the boundary velocities
        return vx,vy
                      
                
    def getCoordinates(self):
        r"""
        Function get all the coordinates of the velocity function space
        :math:`\mathbf{V}`. With the returned coordinates, one cold calculate
        the velocity field in the navier-stokes domain.
        
        The function acquires the DOF coordinates of the vector function space
        by running FEniCS internal function to tabulate all the coordinates.
        This will be in the local coordinate sytem and the will be transformed
        to global coordinate system by **cmGlobal**.
                
        Usage
        -----
        .. code-block:: python
        
            x,y = getCoordinates()
            
        or
        
        .. code-block:: python
        
            xy = getCoordinates()
        
        Parameters
        ----------
        None
        
        Returns
        -------
        xyCoordinates : numpy.ndarray(float64), shape (2,N_DOFs)
                        the global :math:`x,y`-coordinates of the vector 
                        function space  :math:`\mathbf{V}` dofs :math:`N_{DOFs}}`.
            
        Attribute
        ---------
        None changed.
                      
        :First Added:   2013-12-18
        :Last Modified: 2014-03-06
        :Copyright:     Copyright (C) 2014 Lento Manickathan **pHyFlow**
        :Licence:       GNU GPL version 3 or any later version               
        """

        # retrieve the global dof coordinates of the vector function space
        return self.__solver.vectorDOF_coordinates()
        

    def getVorticity(self):
        r"""
        Function to evaluate the vorticity at the probe coordinates defined by
        the probe grid mesh.
        
        Usage
        -----
        .. code-block:: python
        
            wGrid = getVorticity()
            
        Parameters
        ----------
        None        
        
        Returns
        -------
        wGrid : numpy.ndarray(float64), shape (NxProbes,NyProbes)
                the vorticity :math:`\omega` at the  probe grid coordinates
                 in the navier-stokes domain.
                 
        Attributes
        ----------
        __probes : 
                   clear the old probed data and re-probe the vorticity
                   function space.
        
        :First Added:   2013-12-22
        :Last Modified: 2014-02-21
        :Copyright:     Copyright (C) 2014 Lento Manickathan **pHyFlow**
        :Licence:       GNU GPL version 3 or any later version   
        """

        # Should we re-probe the vorticity field
        if self.__reprobeVorticityFlag:
        
            # Remove old vorticity data 
            if self.__probes.number_of_evaluations() > 0:
                self.__probes.clear()
            
            # Probe the data
            self.__probes(self.__solver.vorticity())
            
            # Flow probed, no need to reprobe the current state
            self.__reprobeVorticityFlag = False
        
        
        # Return the data as the shape Nx,Ny of Probes
        return self.__probes.array().reshape(self.__probeGrid['N'][1],self.__probeGrid['N'][0])


    def gTotalInside(self):
        r"""
        Function to calculate the total circulation inside the circulation
        probes coordinates.
        
        Usage
        -----
        .. code-block:: python
        
            gTotalInside = gTotalInside()
            
        Parameters
        ----------
        None        
        
        Returns
        -------
        gTotalInside : float
                       the total circulation inside the probe boundary
                 
        Attributes
        ----------
        __circulationProbes : 
                   clear the old probed data and re-probe the vorticity
                   function space.
        
        :First Added:   2014-03-24
        :Last Modified: 2014-03-24
        :Copyright:     Copyright (C) 2014 Lento Manickathan **pHyFlow**
        :Licence:       GNU GPL version 3 or any later version   
        """

              
        # Remove old vorticity data 
        if self.__circulationProbes.number_of_evaluations() > 0:
            self.__circulationProbes.clear()
            
        # Probe the data
        self.__circulationProbes(self.__solver.u1)

        # Get the data
        vxProbes,vyProbes = self.__circulationProbes.array().T.copy()

        # Make references
        dx = self.__circulationProbes_dl[0]
        dy = self.__circulationProbes_dl[1]
        
        # Integrate the velocity: stokes theorem
        #gTotalInside = _numpy.sum(0.5*(vxProbes[1:]+vxProbes[:-1])*dx + 0.5*(vyProbes[1:]-vyProbes[:-1])*dy)
        gTotalInside = _numpy.sum(vxProbes*dx + vyProbes*dy)
            
        # Return the total circulation inside
        return gTotalInside
        

    def PressureForces(self):
        r"""
        Function to calculate the Pressure gradient forces.
        
        If the free-stream flow is in x-direction (1,0):
            
            Fx : Pressure drag
            Fy : Pressure induced lift
        
        Usage
        -----
        .. code-block:: python
        
            Fp = PressureForces()
            
        Parameters
        ----------
        None
        
        Returns
        -------
        Fp : numpy.ndarray(float64), shape (2,), unit (N)
             the :math:`x,y` component of the pressure forces acting on the
             geometry (defined by the no-slip boundary).
            
        :First Added:   2014-03-05
        :Last Modified: 2014-03-05
        :Copyright:     Copyright (C) 2014 Lento Manickathan **pHyFlow**
        :Licence:       GNU GPL version 3 or any later version     
        """
        return self.__solver.PressureForces()


    def setVelocity(self,vx,vy):
        r"""
        Function to replace the current `(initial)` velocity Field. The input
        requires the :math:`x,y` component of the velocity at the vector
        function space dof coordiantes from .. py:module:: getCoordinates
        
        Usage
        -----
        .. code-block:: python
            
            setVelocity(vx,vy)
        
        
        Parameters
        ----------
        vx : numpy.ndarray(float64), shape (nDOFs,)
             the :math:`x`-component of the velocity at the vector function
             space :math:`\mathbf{V}` DOFs from .. py:module:: getCoordinates.
                     
        vy : numpy.ndarray(float64), shape (nDOFs,)
             the :math:`y`-component of the velocity at the vector function
             space :math:`\mathbf{V}` DOFs from .. py:module:: getCoordinates.
        
        Returns
        -------
        None
        
        Attribute
        ----------
        __solver.u1.vector : numpy.ndarray(float64), shape (2,nDOFs)
                             the :math:`x,y` velocity at the vector function space 
                             :math:`\mathbf{V}` DOFs.
        
        :First Added:   2013-12-18
        :Last Modified: 2014-02-21
        :Copyright:     Copyright (C) 2014 Lento Manickathan **pHyFlow**
        :Licence:       GNU GPL version 3 or any later version        
        """
        
        # vField to assign
        vField = _numpy.vstack((vx,vy)).T.flatten()
        
        # Apply new data
        self.__solver.u1.vector()[:] = vField
                   
               
               
    def __advanceTime(self):
        """
        Function to advance the time.
        
        Usage
        -----
        .. code-block:: python
        
            __advanceTime()
            
        Parameters
        ----------
        None
        
        Returns
        -------
        None returned.
        
        Attribute
        ---------
        __reprobeVorticityFlag
        __solver :
            recalculateVorticityFlag
        __tStep
        __t
              
        :First Added:   2014-02-21
        :Last Modified: 2014-02-21
        :Copyright:     Copyright (C) 2014 Lento Manickathan, **pHyFlow**
        :License:       GNU GPL version 3 or any later version              
        """        
        # advance tStep
        self.__tStep += 1
        
        # advance t
        self.__t = self.__solver.deltaT * self.__tStep
        
        # We need to recalculate the vorticity
        self.__solver.recalculateVorticityFlag = True
        self.__reprobeVorticityFlag = True
       
       
    def __set_variables(self, varName, var):
        """
        Function to check/set the input parameters
        """
        # Check/set geometry data
        if varName == 'geometry':
            if type(var['mesh']) != str:
                raise TypeError("geometry['mesh'] must be of type str. It is %s." % type(var['mesh']))
            
            if type(var['boundaryDomains']) != str:
                raise TypeError("geometry['boundaryDomains'] must be of type str. It is %s." % type(var['boundaryDomains']))
                     
            if type(var['cmGlobal']) != _numpy.ndarray:
                raise TypeError("geometry['cmGlobal'] must be a numpy ndarray. It is %s." % str(type(var['cmGlobal'])))
            if var['cmGlobal'].shape != (2,):
                raise ValueError("%geometry['cmGlobal'] must have shape (2,). It has shape %s." % str(var['cmGlobal'].shape))
            
            if type(var['thetaLocal']) != float and type(var['thetaLocal']) != _numpy.float64:
                raise TypeError("geometry['thetaLocal'] must be a float. It is %s." % str(type(var['thetaLocal'])))
            
            self.__meshFile = var['mesh']
            self.__boundaryDomainsFile = var['boundaryDomains']
            self.__cmGlobal = var['cmGlobal']
            self.__thetaLocal = var['thetaLocal']
            #------------------------------------------------------------------
            
        # Check probegrid parameters
        elif varName == 'probeGrid':
            if type(var['origin']) != _numpy.ndarray:
                raise TypeError("probeGrid['origin'] must be a numpy ndarray. It is %s." % str(type(var['origin'])))
            if var['origin'].shape != (2,):
                raise ValueError("probeGrid['origin']  must have shape (2,). It has shape %s." % str(var['origin'].shape))
            
            if type(var['L']) != _numpy.ndarray:
                raise TypeError("probeGrid['L'] must be a numpy ndarray. It is %s." % str(type(var['L'])))
            if var['L'].shape != (2,):
                raise ValueError("probeGrid['L] must have shape (2,). It has shape %s." % str(var['L'].shape))
            
            if type(var['N']) != _numpy.ndarray:
                raise TypeError("probeGrid['N'] must be a numpy ndarray. It is %s." % str(type(var['N'])))
            if var['N'].shape != (2,):
                raise ValueError("probeGrid['N'] must have shape (2,). It has shape %s." % str(var['N'].shape))

            if type(var['circulationProbes']) != _numpy.ndarray:
                raise TypeError("probeGrid['circulationProbes'] must be a numpy ndarray. It is %s." % str(type(var['circulationProbes'])))
            if var['circulationProbes'].shape[0] != 2:
                raise ValueError("probeGrid['circulationProbes'] must have shape (2,N). It has shape %s." % str(var['circulationProbes'].shape))
          
            self.__probeGrid = var
            #------------------------------------------------------------------

        # Check uMax
        elif varName == 'uMax':
            if type(var) != float and type(var) != _numpy.float64:
                raise TypeError('uMax must be a float. It is %s.' % str(type(var)))
            
            self.__uMax = var
            #------------------------------------------------------------------
        
        elif varName == 'nu':
            if type(var) != float and type(var) != _numpy.float64:
                raise TypeError('nu must be a float. It is %s.' % str(type(var)))
            self.__nu = var
            #------------------------------------------------------------------

        elif varName == 'cfl':
            if type(var) != float and type(var) != _numpy.float64:
                raise TypeError('cfl must be a float. It is %s.' % str(type(var)))
            self.__cfl = var        
            #------------------------------------------------------------------

        elif varName == 'solverParams':
            if var['solver'] not in eulerianOptions.SOLVER['available']:
                raise ValueError(r"solverParams['solver'] must be in [%s]. It is %s" % (eulerianOptions.SOLVER['available'],var['solver'])) 
            self.__solverParams = var
            #------------------------------------------------------------------

        # Set delta T
        elif varName == 'deltaT':
            if type(var) == str:
                if var != 'max':
                    raise TypeError("unknown string parameter: %s. Use 'max'." % var)
                    
            elif type(var) != float and type(var) != _numpy.float64:
                raise TypeError('cfl must be a float. It is %s.' % str(type(var)))
                
            # Set the time-step    
            self.__solver.set_deltaT(var)
            #------------------------------------------------------------------
        
    #--------------------------------------------------------------------------       
    # Attributes

    # Time step size
    @property
    def deltaT(self):
        r"""
        deltaT : float
                 the time step size :math:`\Delta t`.
        """
        return self.__solver.deltaT
        
    @deltaT.setter
    def deltaT(self,deltaTNew):
        self.__set_variables('deltaT',deltaTNew)
                
    # Time step  max
    @simpleGetProperty
    def deltaTMax(self):
        r"""
        deltaTMax : float
                    the maximum allowable time step size :math:`|Delta t_{max}`
        """
        return self.__solver.deltaTMax
                                 
    # cfl number 
    @simpleGetProperty                                 
    def cfl(self):
        r"""
        cfl : float
              the Courant–Friedrichs–Lewy condition stability number,
              :math:`CFL`.
        """
        return self.__cfl

    # cm global coordinates
    @simpleGetProperty
    def cmGlobal(self):
        r"""
        cmGlobal : numpy.ndarray(float64), shape (2,)
                   the :math:`x,y` position of the mesh local reference point
                   (0.,0.) in the global coordinates. 
        """
        return self.__cmGlobal
        
    @simpleGetProperty
    def gTotal(self):
        r"""
        gTotal : float
                 the total circulation :math:`Gamma_t` of the eulerian domain
        """
        return self.__solver.totalCirculation()
                           
    # minimum cell size                           
    @simpleGetProperty
    def hMin(self):
        r"""
        hMin : float
               the minimum mesh cell size 
        """
        return self.__solver.hmin
        
    # viscosity
    @simpleGetProperty                            
    def nu(self):
        """
        nu : float
             the fluid kinematic viscosity :math:`\nu`.
        """
        return self.__nu
                       
    @simpleGetProperty        
    def probeGridMesh(self):
        r"""
        xProbes, yProbes : numpy.ndarray(float64), shape (nxProbes,nyProbes)
                           the *local* :math:`x,y` coordinates of the probe grid mesh.
        """
        # TODO: make sure that this being in local coordinates
        # does not cause any problem
        return self.__xyProbes

    @simpleGetProperty        
    def probeGridParams(self):
        r"""
        probeGridParams : dict
                          dictionary containing all the parameters of the probe
                          grid for extracting the vorticity data.                             
    
                          'origin' : numpy.ndarray(float64), shape (2,)
                                     the :math:`x,y` coordinate of the origin of
                                     the probe mesh in the **local** coordinate
                                     system of the fluid mesh.
                                 
                          'L' : numpy.ndarray(float64), shape (2,)
                                the width :math:`L_x` and the height :math:`L_y`
                                of the probe mesh.
                        
                          'N' : numpy.ndarray(float64), shape (2,)
                                the number of probes :math:`N_x,N_y` in the
                                :math:`x,y` direction.
        """
        return self.__probeGrid
     
    # solver parameters 
    @simpleGetProperty                           
    def solverParams(self):
        r"""
        solverParams : dict
                       the dictionary file containing all the solver parameters.
                       
                       'solver' : str
                                  type of NS solver that should used to solve 
                                  the incompressible laminar problem.
        """
        return self.__solverParams
              
    @simpleGetProperty
    def t(self):
        r"""
        t : float
            the current time of the simulation
        """
        return self.__t                   
                   
    # local theta
    @simpleGetProperty
    def thetaLocal(self):
        """
        thetaLocal : float
                     the local rotational angle :math:`\theta` of the mesh
                     domain. Therefore, the rotation will be done about local                                                                                
                     reference  point (0.,0.), i.e cmGlobal in the global 
                     coordinate system.
        """
        return self.__thetaLocal
        
    @simpleGetProperty
    def tStep(self):
        r"""
        tStep : float
                the current step of the simulation
        """
        return self.__tStep        
                   
    # maximum fluid velocity
    @simpleGetProperty                   
    def uMax(self):
        r"""
        uMax : float
               the maximum fluid velocity :math:`\mathbf{u}_{max}` 
        """
        return self.__uMax
                            


        
    