"""
Navier-Stokes equations for incompressible fluid flow.

Find :math:`\ul{u}`, :math:`p` such that:

.. math::
    \int_{\Omega} \nu\ \nabla \ul{v} : \nabla \ul{u}
    + \int_{\Omega} ((\ul{u} \cdot \nabla) \ul{u}) \cdot \ul{v}
    - \int_{\Omega} p\ \nabla \cdot \ul{v}
    = 0
    \;, \quad \forall \ul{v} \;,

    \int_{\Omega} q\ \nabla \cdot \ul{u}
    = 0
    \;, \quad \forall q \;.
"""

from __future__ import absolute_import
from sfepy import data_dir

filename_mesh ='./pipe_7f_long.mesh'


def poise(ts, coors, bc = None, problem=None):
    x = coors[:,0]
    y = coors[:,1]
    z = coors[:,2]
    return -2*0.459/0.00035**2*(0.00035**2 - (x + 0.0053)**2 - z**2)
        


options = {
    'nls' : 'newton',
    'ls' : 'ls',
    #'post_process_hook' : 'verify_incompressibility',

    # Options for saving higher-order variables.
    # Possible kinds:
    #    'strip' ... just remove extra DOFs (ignores other linearization
    #                options)
    #    'adaptive' ... adaptively refine linear element mesh.
    'linearization' : {
        'kind' : 'strip',
        'min_level' : 1, # Min. refinement level to achieve everywhere.
        'max_level' : 2, # Max. refinement level.
        'eps' : 1e-7, # Relative error tolerance.
    },
}

field_1 = {
    'name' : '3_velocity',
    'dtype' : 'real',
    'shape' : (3,),
    'region' : 'Omega',
    'approx_order' : '1B',
}

field_2 = {
    'name' : 'pressure',
    'dtype' : 'real',
    'shape' : (1,),
    'region' : 'Omega',
    'approx_order' : 1,
}

# Can use logical operations '&' (and), '|' (or).

region_1000 = {
    'name' : 'Omega',
    'select' : 'all',
}

region_0 = {
    'name' : 'Walls',
    'select' : 'vertices of surface -f (r.Outlet +v r.Inlet)', # this works fine for pipe case
# earlier problems were due to the erroneous calculation of flow rate
    'kind' : 'facet',
}
region_1 = {
    'name' : 'Inlet',
    'select' : 'vertices in (y > 0.07999)', #'vertices in (y > 0.07999)',by cinc_0
    'kind' : 'facet',
}
region_2 = {
    'name' : 'Outlet',
    'select' : 'vertices in (y < -0.01999)', #'vertices in (y < -0.01999)',by cinc_1
   'kind' : 'facet',
}



ebc_1 = {
    'name' : 'Inlet',
    'region' : 'Inlet',
#    'dofs' : {'p.0' : 20},
    'dofs' : {'u.1' : -0.459, 'u.[0,2]' : 0.0},
 #-0.459 in case of constant vel at inlet
 # -0.00459 for case 6
}

ebc_3 = {
    'name' : 'Walls',
    'region' : 'Walls',
    'dofs' : {'u.all' : 0.0},
}

#ebc_2 = {
#	'name' : 'Outlet',
#	'region' : 'Outlet',
#	'dofs' : {'p.0' : -20,},
#}

material_1 = {
    'name' : 'fluid',
    'values' : {
        'viscosity' : 0.00365,
        'rho' : 0.001,
    },
}

variable_1 = {
    'name' : 'u',
    'kind' : 'unknown field',
    'field' : '3_velocity',
    'order' : 0,
}
variable_2 = {
    'name' : 'v',
    'kind' : 'test field',
    'field' : '3_velocity',
    'dual' : 'u',
}
variable_3 = {
    'name' : 'p',
    'kind' : 'unknown field',
    'field' : 'pressure',
    'order' : 1,
}
variable_4 = {
    'name' : 'q',
    'kind' : 'test field',
    'field' : 'pressure',
    'dual' : 'p',
}

integral_1 = {
    'name' : 'i1',
    'order' : 2,
}
integral_2 = {
    'name' : 'i2',
    'order' : 3,
}

##
# Stationary Navier-Stokes equations.
equations = {
    'balance' :
    """+ dw_div_grad.i2.Omega( fluid.viscosity, v, u )
       + dw_convect.i2.Omega( v, u )
       - dw_stokes.i1.Omega( v, p ) = 0""",
    #"""+ dw_div_grad.i2.Omega( fluid.viscosity, v, u )
    #   - dw_stokes.i1.Omega( v, p ) = 0""",
    'incompressibility' :
    """dw_stokes.i1.Omega( u, q ) = 0""",
}

solver_0 = {
    'name' : 'ls',
    'kind' : 'ls.scipy_iterative',
    #'kind' : 'ls.petsc' 
    'method' : 'cg',
    'i_max' : 1000,
    'eps_r' : 1e-16,
}

solver_1 = {
    'name' : 'newton',
    'kind' : 'nls.newton',

    'i_max'      : 5, #no change
    'eps_a'      : 1e-8,# original 1e-8, makes a difference. this is the solver stopping criterion
    'eps_r'      : 1.0,
    'macheps'   : 1e-16,
    'lin_red'    : 1e-6, # Linear system error < (eps_a * lin_red).# no change
    'ls_red'     : 0.1,
    'ls_red_warp' : 0.001,# changed from 0.001 no change
    'ls_on'      : 0.99999,# original 0.99999 no change
    'ls_min'     : 1e-8, # changed form 1e-5, 1e-3 no change, 1e-8 no change
    'check'     : 0,
    'delta'     : 1e-8, # no change
}

functions = {
    'poise' : (poise,),
}

def verify_incompressibility( out, problem, state, extend = False ):
    """This hook is normally used for post-processing (additional results can
    be inserted into `out` dictionary), but here we just verify the weak
    incompressibility condition."""
    from sfepy.base.base import nm, output, assert_

    vv = problem.get_variables()
    one = nm.ones( (vv['p'].field.n_nod,), dtype = nm.float64 )
    vv['p'].set_data(one)
    zero = problem.evaluate('dw_stokes.i1.Omega( u, p )', p=one, u=vv['u']())
    output('div( u ) = %.3e' % zero)

    assert_(abs(zero) < 1e-14)

    return out


import utils as utils

cinc = getattr(utils, 'cinc_cylinder')

functions = {
    'cinc_0' : (lambda coors, domain=None: cinc(coors, 0),),
    'cinc_1' : (lambda coors, domain=None: cinc(coors, 1),),
}