# Path: dolfin_navier_scipy/residual_checks.py
# def get_imex_res(V=None, outflowds=None, gradvsymmtrc=True, nu=None,
#                  implscheme='crni', explscheme='abtw'):
#     """ define the residual for an IMEX/AB2 time discretization
# 
#     """
#     if not implscheme == 'crni':
#         raise NotImplementedError()
# 
#     if explscheme == 'abtw':
#         def convform(cvo=None, cvt=None, phi=None):
#             return (1.5*inner(dolfin.dot(cvo, nabla_grad(cvo)), phi)*dx -
#                     .5*inner(dolfin.dot(cvt, nabla_grad(cvt)), phi)*dx)
# 
#     elif explscheme == 'heun':
#         def convform(cvo=None, cvt=None, phi=None):
#             return (.5*inner(dolfin.dot(cvo, nabla_grad(cvo)), phi)*dx +
#                     .5*inner(dolfin.dot(cvt, nabla_grad(cvt)), phi)*dx)
# 
#     elif explscheme == 'eule':
#         def convform(cvo=None, cvt=None, phi=None):
#             return inner(dolfin.dot(cvo, nabla_grad(cvo)), phi)*dx
# 
#     def imex_res(vel, pres, dt, lastvel=None, othervel=None, phi=None):
#         if phi is None:
#             phi = dolfin.TestFunction(V)
# 
#         diffvel = .5*(vel+lastvel)  # Crank-Nicolson
#         diffrm = nu*inner(grad(diffvel)+grad(diffvel).T, grad(phi))*dx
#         if gradvsymmtrc:
#             nvec = dolfin.FacetNormal(V.mesh())
#             diffrm = diffrm - (nu*inner(grad(diffvel).T*nvec, phi))*outflowds
#         cnvfrm = convform(cvo=lastvel, cvt=othervel, phi=phi)
# 
#         pfrm = -inner(pres, div(phi))*dx
#         dtprt = 1./dt*dolfin.assemble(inner(vel, phi)*dx) \
#             - 1./dt*dolfin.assemble(inner(lastvel, phi)*dx)
#         res = dolfin.assemble(diffrm+cnvfrm+pfrm) + dtprt
# 
#         # import numpy as np
#         # nfc_c = dolfin.assemble(cnvfrm).get_local()
#         # print(np.linalg.norm(nfc_c), nfc_c[0], nfc_c.size)
#         # print('debggng')
#         return res
# 
#     return imex_res


import numpy as np
import matplotlib.pyplot as plt
import json
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_imex_res

    # function of ones at the lift/drag boundary
    phionevec = np.zeros((femp['V'].dim(), 1))
    phionevec[femp['mvwbcinds'], :] = 1.
    phione = dolfin.Function(femp['V'])
    phione.vector().set_local(phionevec)
    pickx = dolfin.as_matrix([[1., 0.], [0., 0.]])
    picky = dolfin.as_matrix([[0., 0.], [0., 1.]])
    pox = pickx*phione
    poy = picky*phione

    # function of the tangential vector at the lift/drag boundary
    phitwovec = np.zeros((femp['V'].dim(), 1))
    phitwovec[femp['mvwbcinds'], 0] = femp['mvwbcvals']
    phitwo = dolfin.Function(femp['V'])
    phitwo.vector().set_local(phitwovec)

    # getld = dnsps.LiftDragSurfForce(V=femp['V'], nu=nu,
    #                                 phione=phione, phitwo=phitwo,
    #                                 outflowds=femp['outflowds'],
    #                                 ldds=femp['liftdragds'])

    # L = femp['charlen']  # characteristic length = 2*Radius

    a_1 = dolfin.Point(0.15, 0.2)
    a_2 = dolfin.Point(0.25, 0.2)

    reschkdict = dict(V=femp['V'], gradvsymmtrc=True,
                      outflowds=femp['outflowds'], nu=nu)
    euleres = get_imex_res(explscheme='eule', **reschkdict)

# Path: dolfin_navier_scipy/residual_checks.py
# def get_imex_res(V=None, outflowds=None, gradvsymmtrc=True, nu=None,
#                  implscheme='crni', explscheme='abtw'):
#     """ define the residual for an IMEX/AB2 time discretization
# 
#     """
#     if not implscheme == 'crni':
#         raise NotImplementedError()
# 
#     if explscheme == 'abtw':
#         def convform(cvo=None, cvt=None, phi=None):
#             return (1.5*inner(dolfin.dot(cvo, nabla_grad(cvo)), phi)*dx -
#                     .5*inner(dolfin.dot(cvt, nabla_grad(cvt)), phi)*dx)
# 
#     elif explscheme == 'heun':
#         def convform(cvo=None, cvt=None, phi=None):
#             return (.5*inner(dolfin.dot(cvo, nabla_grad(cvo)), phi)*dx +
#                     .5*inner(dolfin.dot(cvt, nabla_grad(cvt)), phi)*dx)
# 
#     elif explscheme == 'eule':
#         def convform(cvo=None, cvt=None, phi=None):
#             return inner(dolfin.dot(cvo, nabla_grad(cvo)), phi)*dx
# 
#     def imex_res(vel, pres, dt, lastvel=None, othervel=None, phi=None):
#         if phi is None:
#             phi = dolfin.TestFunction(V)
# 
#         diffvel = .5*(vel+lastvel)  # Crank-Nicolson
#         diffrm = nu*inner(grad(diffvel)+grad(diffvel).T, grad(phi))*dx
#         if gradvsymmtrc:
#             nvec = dolfin.FacetNormal(V.mesh())
#             diffrm = diffrm - (nu*inner(grad(diffvel).T*nvec, phi))*outflowds
#         cnvfrm = convform(cvo=lastvel, cvt=othervel, phi=phi)
# 
#         pfrm = -inner(pres, div(phi))*dx
#         dtprt = 1./dt*dolfin.assemble(inner(vel, phi)*dx) \
#             - 1./dt*dolfin.assemble(inner(lastvel, phi)*dx)
#         res = dolfin.assemble(diffrm+cnvfrm+pfrm) + dtprt
# 
#         # import numpy as np
#         # nfc_c = dolfin.assemble(cnvfrm).get_local()
#         # print(np.linalg.norm(nfc_c), nfc_c[0], nfc_c.size)
#         # print('debggng')
#         return res
# 
#     return imex_res


# Path: tests/tdp_2D_simu.py
import numpy as np
import matplotlib.pyplot as plt
import json
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_imex_res


    # function of ones at the lift/drag boundary
    phionevec = np.zeros((femp['V'].dim(), 1))
    phionevec[femp['mvwbcinds'], :] = 1.
    phione = dolfin.Function(femp['V'])
    phione.vector().set_local(phionevec)
    pickx = dolfin.as_matrix([[1., 0.], [0., 0.]])
    picky = dolfin.as_matrix([[0., 0.], [0., 1.]])
    pox = pickx*phione
    poy = picky*phione

    # function of the tangential vector at the lift/drag boundary
    phitwovec = np.zeros((femp['V'].dim(), 1))
    phitwovec[femp['mvwbcinds'], 0] = femp['mvwbcvals']
    phitwo = dolfin.Function(femp['V'])
    phitwo.vector().set_local(phitwovec)

    # getld = dnsps.LiftDragSurfForce(V=femp['V'], nu=nu,
    #                                 phione=phione, phitwo=phitwo,
    #                                 outflowds=femp['outflowds'],
    #                                 ldds=femp['liftdragds'])

    # L = femp['charlen']  # characteristic length = 2*Radius

    a_1 = dolfin.Point(0.15, 0.2)
    a_2 = dolfin.Point(0.25, 0.2)

    reschkdict = dict(V=femp['V'], gradvsymmtrc=True,
                      outflowds=femp['outflowds'], nu=nu)
    euleres = get_imex_res(explscheme='eule', **reschkdict)

    heunres = get_imex_res(explscheme='heun', **reschkdict)