Local solver gives different results for two equivalent definitions of "ds"-measure

[Tormod Landet]

The following ds-measure should be equivalent to the normal ds measure:

marker = FacetFunction("size_t", mesh)
ds0 = Measure("ds", domain=mesh, subdomain_data=marker, subdomain_id=0)

However, if I use "ds" for the bilinear form and "ds0" for the linear form (or vice versa) I get different results than if I use either of them for both forms

from dolfin import *

mesh = UnitTriangleMesh()
marker = FacetFunction("size_t", mesh)
ds0 = Measure("ds", domain=mesh, subdomain_data=marker, subdomain_id=0)

Vx = VectorFunctionSpace(mesh, 'DG', 1)
Vs = FunctionSpace(mesh, 'DGT', 1)

u = TrialFunction(Vx)
v = TestFunction(Vs)
n = FacetNormal(mesh)
a = dot(u, n)*v*ds
L = v*ds0

uh1 = Function(Vx)    
ls = LocalSolver(a, L)
ls.solve_local_rhs(uh1)
print norm(uh1)

uh2 = Function(Vx)
A, b = assemble_system(a, L)
solve(A, uh2.vector(), b)
print norm(uh2)

assert errornorm(uh1, uh2, degree_rise=0) < 1e-10

Changing "a = dot(u, n)vds" to be "a = dot(u, n)vds0" or using "ds" in both terms makes the code work as expected. ```

[Tormod Landet]

Just to confirm that the error must be in the LocalSolver and not in the LocalAssembler, here is a comparison with a quick and dirty Python implementation of the local solver:

from dolfin import *

mesh = UnitTriangleMesh()
marker = FacetFunction("size_t", mesh)
ds0 = Measure("ds", domain=mesh, subdomain_data=marker, subdomain_id=0)

Vx = VectorFunctionSpace(mesh, 'DG', 1)
Vs = FunctionSpace(mesh, 'DGT', 1)

u = TrialFunction(Vx)
v = TestFunction(Vs)
n = FacetNormal(mesh)
a = dot(u, n)*v*ds
L = v*ds0

class PythonLocalSolver(object):
    def __init__(self, a, L):
        self.a = a
        self.L = L

    def solve_local_rhs(self, f):
        from numpy.linalg import solve as npsolve
        V = f.function_space()
        mesh = V.mesh()        

        for cell in cells(mesh):
            Ai = assemble_local(self.a, cell)
            bi = assemble_local(self.L, cell)
            ui = npsolve(Ai, bi)
            dofs = V.dofmap().cell_dofs(cell.index())
            for i, dof in enumerate(dofs):
                f.vector()[dof] = ui[i]

# Dolfin local solver
uhD = Function(Vx)    
ls = LocalSolver(a, L)
ls.solve_local_rhs(uhD)
print 'uhD:', uhD.vector().array()

# Python local solver
uhP = Function(Vx)    
ls = PythonLocalSolver(a, L)
ls.solve_local_rhs(uhP)
print 'uhP:', uhP.vector().array()

assert errornorm(uhD, uhP, degree_rise=0) < 1e-10

The pure python version gives the expected result no matter which of the equivalent "ds" measures are used in the linear/bilinear form

[Tormod Landet]

Looking at LocalSolver.cpp makes the error quite apparent.

// Extract cell_domains etc from left-hand side form

This needs to be fixed to use the cell domains from the right-hand-side with the linear form

[Tormod Landet]

This is fixed by https://bitbucket.org/fenics-project/dolfin/pull-requests/281

[Anders Logg (Chalmers)]

Removing milestone: 1.7 (automated comment)

[Prof Garth Wells]

Fix issue 696 - Local solver disregards RHS form domains

→ <<cset 76a27a559238>>

[Prof Garth Wells]

Add test to see that we fixed issue 696

We create a facet function and a ds Measure using this facet function. We need not mark any facets and only use marker 0 (which is the default). This was enough to make the previous implementation assemble the form incorrectly. With the fixed LocalSolver RHS domains are handled correctly and the correct local tensor is assembled leading to the test not failing.

→ <<cset b195286ab85f>>