Counter-integrating does not lead to same results

variational formulation
,
1

Hello,

I want to compute

\int_{\Omega} b \nabla u \, v \,\mathrm{d}x = \int_{\Omega} f \, v \, \mathrm{d}x,

where \Omega = [0,1]^2 with inflow boundary conditions u_{|\partial \Omega^{-}} = 1 for \partial \Omega^{-} = \{0\} \times [0,1] and \partial \Omega^{+} = \partial \Omega \setminus \partial \Omega^{-}.

By Greene’s identity

\int_{\Omega} b \nabla u \, v \,\mathrm{d}x = - \int_{\Omega} u \, \mathrm{div} b \, v \, \mathrm{d}x - \int_{\Omega} u \, b \nabla v \,\mathrm{d}x + \int_{\partial \Omega^{+}} b \cdot n \, u v \, \mathrm{d}s,

having used that the integral vanishes on \partial \Omega^{-}. Using either the left or right-hand side to numerically solve the problem should give the same result, but that is not the case in my code. Can someone point out what I am doing wrong?

from dolfin import *
import numpy as np

nx = ny = 20 ; degree = 1

mesh = UnitSquareMesh(nx,ny)
n = FacetNormal(mesh)
V = FunctionSpace(mesh, "CG", degree)
Vf = VectorFunctionSpace(mesh, "CG", degree=degree+1, dim=2)

f = Expression("cos(x[0] * M_PI) + sin(x[1] * M_PI)", degree=degree+1)
beta = Expression(('cos(x[0])', 'sin(2*x[1] - .2)'), degree=degree+1)
betaInterp = interpolate(beta, Vf)

boundaries = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
boundaries.set_all(0)

class Outflow(SubDomain):
    def inside(self, x, on_boundary):
        tol = 1e-13
        if x[0] > tol and on_boundary:
            return True

class Inflow(SubDomain):
    def inside(self, x, on_boundary):
        tol = 1e-13
        if near(x[0], 0, tol) and on_boundary:
            return True

outflow = Outflow()
outflow.mark(boundaries, 1)

inflow = Inflow()
inflow.mark(boundaries, 2)

ds = Measure("ds", subdomain_data=boundaries)

u = TrialFunction(V) ; v = TestFunction(V)

a_v1 = dot(betaInterp, grad(u)) * v * dx

a_v2 = - u * dot(betaInterp, grad(v)) * dx \
    - div(betaInterp) * u *  v * dx \
    + dot(betaInterp, n) * u * v * ds(1)

f_rhs = f * v * dx

A_v1 = assemble(a_v1) ; A_v2 = assemble(a_v2) ; F = assemble(f_rhs)

bc = DirichletBC(V, 1.0, inflow)
bc.apply(A_v1, F)
bc.apply(A_v2)

u1 = Function(V) ; solve(A_v1, u1.vector(), F)
u2 = Function(V) ; solve(A_v2, u2.vector(), F)

print(errornorm(u1,u2)) #output is 18.217744114223798
print(np.where( np.abs(A_v1.array() - A_v2.array()) > 1e-14 )) # output is (array([  2,   2, 231, 231]), array([  0,   2, 210, 231]))

The last line shows where the two assembled matrices differ, which are only 4 entries.

Apologies if this question is basic and/or has already been asked in a different form, I couldn’t find any answer to it.
Thanks in advance.

2

After some fiddling, setting the boundary conditions as follows

class Inflow(dol.SubDomain):
    def inside(self, x, on_boundary):
        if dol.near(x[0],0):
            return True
        return False

class Outflow(dol.SubDomain):
    def inside(self, x, on_boundary):
        if dol.near(x[0], 1.) or dol.near(x[1], 0.) or dol.near(x[1],1.):
            return True
        return False

solves the inconsistency problem.