- edited description
Local solver gives different results for two equivalent definitions of "ds"-measure
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. ```
Comments (7)
-
reporter -
reporter 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
-
reporter 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
-
reporter This is fixed by https://bitbucket.org/fenics-project/dolfin/pull-requests/281
-
- removed milestone
Removing milestone: 1.7 (automated comment)
-
- changed status to resolved
Fix issue 696 - Local solver disregards RHS form domains
→ <<cset 76a27a559238>>
-
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>>
- Log in to comment