# Path: dolfin_navier_scipy/residual_checks.py
# def get_steady_state_res(V=None, outflowds=None, gradvsymmtrc=True, nu=None):
# 
#     def steady_state_res(vel, pres, phi=None):
#         if phi is None:
#             phi = dolfin.TestFunction(V)
# 
#         cnvfrm = inner(dolfin.dot(vel, nabla_grad(vel)), phi)*dx
#         diffrm = nu*inner(grad(vel)+grad(vel).T, grad(phi))*dx
#         if gradvsymmtrc:
#             nvec = dolfin.FacetNormal(V.mesh())
#             diffrm = diffrm - (nu*inner(grad(vel).T*nvec, phi))*outflowds
# 
#         pfrm = inner(pres, div(phi))*dx
#         res = dolfin.assemble(diffrm+cnvfrm-pfrm)
#         return res
# 
#     return steady_state_res


import numpy as np
import dolfin
import dolfin_navier_scipy.stokes_navier_utils as snu
import dolfin_navier_scipy.dolfin_to_sparrays as dts
import dolfin_navier_scipy.problem_setups as dnsps
from dolfin_navier_scipy.residual_checks import get_steady_state_res
    ddir = 'data/'
    data_prfx = problem+'{2}_mesh{0}_Re{1}'.format(meshlvl, femp['Re'], scheme)

    # ## Parameters for the benchmark values
    Um = charvel  # (we alread scale the inflow parabola accordingly)
    L = femp['charlen']  # characteristic length
    NP, NV = stokesmatsc['J'].shape
    print('NV + NP : {0} + {1} = {2}'.format(NV, NP, NV+NP))

    soldict = stokesmatsc  # containing A, J, JT
    # soldict.update(femp)  # adding V, Q, invinds, diribcs
    soldict.update(invinds=femp['invinds'], V=femp['V'], Q=femp['Q'])
    soldict.update(fv=rhsd['fv'], fp=rhsd['fp'],
                   N=meshlvl, nu=nu,
                   verbose=True,
                   return_vp=True,
                   get_datastring=None,
                   dbcinds=femp['dbcinds'], dbcvals=femp['dbcvals'],
                   data_prfx=ddir+data_prfx,
                   paraviewoutput=ParaviewOutput,
                   vfileprfx=proutdir+'vel_',
                   pfileprfx=proutdir+'p_')

#
# compute the uncontrolled steady state Navier-Stokes solution
#
    vp_ss_nse = snu.solve_steadystate_nse(**soldict)
    vss, dynpss = dts.expand_vp_dolfunc(vc=vp_ss_nse[0], pc=vp_ss_nse[1],
                                        **femp)
    steady_state_res = \

# Path: dolfin_navier_scipy/residual_checks.py
# def get_steady_state_res(V=None, outflowds=None, gradvsymmtrc=True, nu=None):
# 
#     def steady_state_res(vel, pres, phi=None):
#         if phi is None:
#             phi = dolfin.TestFunction(V)
# 
#         cnvfrm = inner(dolfin.dot(vel, nabla_grad(vel)), phi)*dx
#         diffrm = nu*inner(grad(vel)+grad(vel).T, grad(phi))*dx
#         if gradvsymmtrc:
#             nvec = dolfin.FacetNormal(V.mesh())
#             diffrm = diffrm - (nu*inner(grad(vel).T*nvec, phi))*outflowds
# 
#         pfrm = inner(pres, div(phi))*dx
#         res = dolfin.assemble(diffrm+cnvfrm-pfrm)
#         return res
# 
#     return steady_state_res


# Path: tests/steadystate_schaefer-turek_2D-1.py
import numpy as np
import dolfin
import dolfin_navier_scipy.stokes_navier_utils as snu
import dolfin_navier_scipy.dolfin_to_sparrays as dts
import dolfin_navier_scipy.problem_setups as dnsps
from dolfin_navier_scipy.residual_checks import get_steady_state_res

    ddir = 'data/'
    data_prfx = problem+'{2}_mesh{0}_Re{1}'.format(meshlvl, femp['Re'], scheme)

    # ## Parameters for the benchmark values
    Um = charvel  # (we alread scale the inflow parabola accordingly)
    L = femp['charlen']  # characteristic length
    NP, NV = stokesmatsc['J'].shape
    print('NV + NP : {0} + {1} = {2}'.format(NV, NP, NV+NP))

    soldict = stokesmatsc  # containing A, J, JT
    # soldict.update(femp)  # adding V, Q, invinds, diribcs
    soldict.update(invinds=femp['invinds'], V=femp['V'], Q=femp['Q'])
    soldict.update(fv=rhsd['fv'], fp=rhsd['fp'],
                   N=meshlvl, nu=nu,
                   verbose=True,
                   return_vp=True,
                   get_datastring=None,
                   dbcinds=femp['dbcinds'], dbcvals=femp['dbcvals'],
                   data_prfx=ddir+data_prfx,
                   paraviewoutput=ParaviewOutput,
                   vfileprfx=proutdir+'vel_',
                   pfileprfx=proutdir+'p_')

#
# compute the uncontrolled steady state Navier-Stokes solution
#
    vp_ss_nse = snu.solve_steadystate_nse(**soldict)
    vss, dynpss = dts.expand_vp_dolfunc(vc=vp_ss_nse[0], pc=vp_ss_nse[1],
                                        **femp)
    steady_state_res = \

        get_steady_state_res(V=femp['V'], gradvsymmtrc=True,