Incorrect integration for exterior_facet with Nedelec elements
Affects | Status | Importance | Assigned to | Milestone | |
---|---|---|---|---|---|
DOLFIN |
Invalid
|
Undecided
|
Unassigned |
Bug Description
Hi,
I expect the following code to be equivalent to a line integral over the edge of the square domain. However it seems as if integration is actually happening in the cells that touch the edge of the domain. (This can be seen by the values for x and y printed out during the assembly)
Do I understand this correctly or is my implementation wrong or is this a bug?
import dolfin as D
mesh = D.UnitSquare(2,2)
V = D.FunctionSpace
bd = D.DomainBoundary()
# mf = D.FacetFunction
mf = D.MeshFunction(
mf.set_all(0)
bd.mark(mf, 2)
dss = D.ds[mf]
print mf.array()
code = """
class MyFunc : public Expression
{
public:
boost:
MyFunc() : Expression(2), mesh(static_
void eval(Array<double>& values, const Array<double>& x, const ufc::cell& ufc_cell) const
{
Cell cell(*mesh, ufc_cell.index);
Point n = cell.normal(
printf ("x: %f, y: %f, n.x: %f, n.y: %f\\n", x[0], x[1], n.x(), n.y());
values[0] = x[0];
values[1] = x[1];
}
};
"""
ex = D.Expression(code, element=
ex.mesh = mesh
v = D.TestFunction(V)
L = D.dot(v, ex)*dss(2)
b = D.assemble(L)
print b.array() # expect this to only have values for edge basis functions defined of the edges
Regards,
Renier Marchand
The printouts from your expression, show that the expression is evaluated at the midpoints of each edge for each cell, e.g. these lines for one particular cell:
x: 1.000000, y: 0.250000, n.x: 0.000000, n.y: -1.000000
x: 0.750000, y: 0.250000, n.x: 0.000000, n.y: -1.000000
x: 0.750000, y: 0.000000, n.x: 0.000000, n.y: -1.000000
What happens is that:
1) Your expression is interpolated into the local fem basis on a cell, which involves the above evaluations on all edges of the cell.
2) Integration over the facet is performed, using the dofs on the cell computed in step 1.
So the prints look fine, and I'm closing this bug as invalid. If you get wrong results in the integration, you need to justify that otherwise and reopen the bug.