wrong interpolation for enriched element spaces
Interpolation to enriched spaces doesn't work as expected....
from dolfin import *
mesh = UnitSquareMesh(2, 2)
V = FunctionSpace(mesh, "CG", 2)
B = FunctionSpace(mesh, "Bubble", 3)
U=V+B
e=Expression("x[0]*x[1]")
u=interpolate(e,U); plot(u) #looks fine
W=FunctionSpace(refine(refine(refine(mesh))),'CG',1)
f=interpolate(e,W); plot(f) #looks fine
f=interpolate(u,W); plot(f) #not good
interactive()
Comments (18)
-
-
reporter Forget the plots....
There is something wrong - even with
e=Expression("x[0]")
which is exactly represented in all the spaces.print u.vector().max(), f.vector().max()
gives1.0 1.666015625
That is not correct for nodal DOFs interpolation.(assuming the Bubble function has nodal DOF as value at the cell center for example)
-
- changed status to invalid
-
- changed status to open
Ok, I see your point. I'll check this.
-
reporter the
f
was just a way to see the bumps inu
....I will modify the text - the problem is already the Expression interpolation to the enriched space.
-
Following code shows that Bubble DOFs are twice as expected
from dolfin import * mesh = UnitSquareMesh(2, 2) V = FunctionSpace(mesh, 'Lagrange', 2) B = FunctionSpace(mesh, 'Bubble', 3) U = V + B e = Expression('x[0]*x[1]', domain=mesh) u = interpolate(e, U) gdim = mesh.geometry().dim() print "P2 dofs" coords = V.dofmap().tabulate_all_coordinates(mesh) gdim = mesh.geometry().dim() for i in range(coords.size / gdim): print coords[2*i:2*i+gdim], e(coords[2*i:2*i+gdim]), u(coords[2*i:2*i+gdim]) print "B3 dofs" coords = B.dofmap().tabulate_all_coordinates(mesh) gdim = mesh.geometry().dim() for i in range(coords.size / gdim): print coords[2*i:2*i+gdim], e(coords[2*i:2*i+gdim]), u(coords[2*i:2*i+gdim])
-
- changed component to dofmaps
- changed version to dev
- marked as critical
- changed milestone to 1.6
-
assigned issue to
-
reporter - edited description
-
Enriched elements are not nodal because the tabulation happens by splitting together the tabulation of the sub elements. To make this work, the basis for the enriched element needs to be reorthogonalised I think.
-
reporter Ok - thanks for the explanation.
1) if the
interpolate
function does straight forward nodal dof interpolationand
2) if the enriched element basis is just the (nodal) basis of those two subspaces put together (=> i.e. sum of all basis functions of the enriched space is > 1 inside each cell)
then the result "makes sense"....
-
reporter some more observations:
wrong interpolation to the enriched space (now expected since basis of
V+B
is not orthogonal with respect to nodal evaluation)mesh = UnitSquareMesh(1, 1, 'crossed') V = FunctionSpace(mesh, "CG", 2) B = FunctionSpace(mesh, "Bubble", 3) U=V+B e=Expression("x[0]*x[1]", domain=mesh) p=Point(0.75,0.5) print e(p) # == 0.375 the correct value Pe=project(e,U); print Pe(p) # == 0.375000016333 Ie=interpolate(e,U); print Ie(p) # == 0.7265625 expected wrong value
assuming the projection
Pe
to the enriched space is fine, if I want to see the function inCG1
on refined meshW=FunctionSpace(refine(mesh),'CG',1) f=project(Pe,W); print f(p) # == 0.541071444256 wrong f=interpolate(Pe,W); print f(p) # == 0.375000016333
but now the projection from the enriched space to CG1 is wrong, while interpolation from the enriched space is fine...
-
@jxh
but now the projection from the enriched space to CG1 is wrong, while interpolation from the enriched space is fine...
I'd guess you're just lucky with the interpolation because
p
is on facet midpoint andf(p)
is just average between two CG1 nodal values:0.5*(Pe(0.75, 0.25)+Pe(0.75, 0.75)) == 0.375
.On the other hand, I wouldn't expect the projection to preserve value of
f(p)
. It isL^2
projection of 3rd-order polynomial to Lagrange degree 1. Please, correct me if I'm wrong. -
I checked code generated for
element = EnrichedElement(FiniteElement('Lagrange', Domain(Cell('interval', 1)), 2, None), FiniteElement('Bubble', Domain(Cell('interval', 1)), 3, None))
by inspecting generated code and
#include <vector> #include <iostream> #include "form.h" int main() { form_finite_element_0 element; form_dofmap_0 dofmap; std::size_t gdim = element.geometric_dimension(); std::size_t dim = element.space_dimension(); std::vector<double> values(dim); std::vector<double> vertex_coords{0.0, 1.0}; std::vector<double> dof_coords(gdim*dim); std::vector<double> x(gdim); dofmap.tabulate_coordinates(dof_coords.data(), vertex_coords.data()); for (std::size_t i = 0; i < dim; ++i) { x.assign(&dof_coords[gdim*i], &dof_coords[gdim*(i+1)]); element.evaluate_basis_all(values.data(), x.data(), vertex_coords.data(), 0); std::cout << "x: "; for (auto xi : x) std::cout << xi << "; "; std::cout << "values: "; for (auto v : values) std::cout << v << "; "; std::cout << std::endl; } return 0; }
and everything seems correct there. This is just
ufc::finite_elment
andufc::dofmap
. I haven't looked yet intoufc::form
but I expect that everything is correct on form compiler side and assembly.Summing up: there's a wrong assumption in DOLFIN interpolation code that DOF functionals are conjugate to (possibly enriched) basis, i.e.
L_i (\psi_j) = \delta_{ij}
forL_i
DOF functional,\psi_j
basis function. Note that there's no notion of orthogonality there. As a solution I'd propose just disabling interpolation where this relation is not satisfied. Question is how to detect this. Later we could consider recomputingL_i
or\psi_j
to match, maybe on FFC side. Any thoughts, @martinal, @meg? -
reporter @blechta
I'd guess you're just lucky with the interpolation because p is on facet midpoint and f(p) is just average between two CG1 nodal values: 0.5*(Pe(0.75, 0.25)+Pe(0.75, 0.75)) == 0.375.
On the other hand, I wouldn't expect the projection to preserve value of f(p). It is L^2 projection of 3rd-order polynomial to Lagrange degree 1. Please, correct me if I'm wrong.
To make it simple - with constant expression at point (0.5,0.5):
from dolfin import * mesh = UnitSquareMesh(1, 1, 'crossed') V = FunctionSpace(mesh, "CG", 2) B = FunctionSpace(mesh, "Bubble", 3) U=V+B e=Expression("1.0", domain=mesh) p=Point(0.5,0.5) print e(p) # 1.0 Pe=project(e,U); print Pe(p) # 0.99999998307 Ie=interpolate(e,U); print Ie(p) # 1.0 W=FunctionSpace(refine(mesh),'CG',1) f=project(Pe,W); print f(p) # 1.44999976228 f=interpolate(Pe,W); print f(p) # 0.99999998307
other observations:
- the problem remains if the space V is CG1
- if the mesh is not refined for the subsequent projection then the result is ok
-
@jxh's example works with
W=FunctionSpace(mesh,'CG',1)
but fails withW=FunctionSpace(Mesh(mesh),'CG',1)
. The first one collects the DOFs directly from vector when assembling rhs linear form while the latter one callsrestrict_as_ufc_function(w, element, dolfin_cell, vertex_coordinates, ufc_cell);
which must hence be broken also. -
"Fix" (disabling faulty features) in https://bitbucket.org/fenics-project/ffc/pull-request/19. This effectively disables interpolation to enriched space, projection from enriched to another mesh and using expressions and constants in Dirichlet BCs on enriched space. The rest (i.e. projection to enriched, interpolation from enriched to non-enriched on the same mesh and Functions in Dirichlet values on enriched space) should be working correctly.
-
- changed status to resolved
Incorrect code is not generated now. For correct reimplementation watch https://bitbucket.org/fenics-project/ffc/issue/69.
-
- removed milestone
Removing milestone: 1.6 (automated comment)
- Log in to comment
I'd say this is expected. Note that
plot(u)
is just interpolation to vertices of functionu
which is nodal (at P2 + B3 DOFs; note that all the DOFs are nodal in DOLFIN) interpolation ofe
.Thus
f=interpolate(u,W); plot(f)
plotsu
much accurately thanplot(u)
. This is known behaviour in DOLFIN mainly caused by restrictions given by VTK.