assemble on mixed meshview-forms returns wrong results

Issue #1123 new
Johannes Haubner created an issue

We found a bug when assembling forms that contain functions/measures from different meshes using meshview. When assembling this form with assemble, the resulting vector are incorrect. If one uses assemble_mixed the results are correct. A minimal example of the problem is shown below. We see two possible solutions to this issue:

  1. Throw an error message when calling assemble with forms that contain different mesh-view meshes.
  2. Add support in assemble for mixed forms.

We believe that option 2 is preferable over 1, since then helper functions, e.g. solve should automatically work also mixed mesh-view problems.

Minimal example:

from dolfin import *
from matplotlib import pyplot

mesh = UnitSquareMesh(24, 24)

class DesignBoundary(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and abs(x[0]) < DOLFIN_EPS

# Create mesh view
marker = MeshFunction("size_t", mesh, 1)
marker.set_all(0)
design_boundary = DesignBoundary()
design_boundary.mark(marker, 1)
design_mesh = MeshView.create(marker, 1)

# Create custom measure
dxb = Measure("dx", domain = design_mesh)

# initialize function spaces
W = FunctionSpace(mesh, "CG", 1)
B = FunctionSpace(design_mesh, "CG", 1)

# define bilinear form
b = project(Constant(1.0), B)
w = TrialFunction(W)
psiw = TestFunction(W)
a = w*psiw*dx
L = b*psiw*dxb
A = assemble(a)
lhs = assemble_mixed(L)  # yields correct results
lhs2 = assemble(L)   # yields wrong results
w_ = Function(W)
solve(A, w_.vector(), lhs2)

plot(w_)
pyplot.savefig("w_meshview.png")
pyplot.cla()

We demonstrate the problem in this example, where we project a function b that lives on a part of the boundary of the mesh to a function w on the whole mesh. Using assemble for the right hand side of the pde yields a wrong result (first picture), whereas mixed_assemble yields the correct result (second picture).

We tried to fix it by replacing the Assembler by the MixedAssembler ( https://bitbucket.org/JohannesHaubner/dolfin/commits/5a2ae044359971dfb3b88be4aeae50b400bc7148 ). It worked for our minimal example and also for some other examples that we tried.

Comments (3)

  1. Log in to comment