wrong interpolation for enriched element spaces

[Jaroslav Hron]

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()

[Jan Blechta]

I'd say this is expected. Note that plot(u) is just interpolation to vertices of function u which is nodal (at P2 + B3 DOFs; note that all the DOFs are nodal in DOLFIN) interpolation of e.

Thus f=interpolate(u,W); plot(f) plots u much accurately than plot(u). This is known behaviour in DOLFIN mainly caused by restrictions given by VTK.

[Jaroslav Hron]

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() gives 1.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)

[Jan Blechta]

Ok, I see your point. I'll check this.

[Jaroslav Hron]

the f was just a way to see the bumps in u....

I will modify the text - the problem is already the Expression interpolation to the enriched space.

[Jan Blechta]

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])

[Lawrence Mitchell]

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.

[Jaroslav Hron]

Ok - thanks for the explanation.

1) if the interpolate function does straight forward nodal dof interpolation

and

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"....

[Jaroslav Hron]

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 in CG1 on refined mesh

W=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...

[Jan Blechta]

@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 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.

[Jan Blechta]

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 and ufc::dofmap. I haven't looked yet into ufc::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} for L_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 recomputing L_i or \psi_j to match, maybe on FFC side. Any thoughts, @martinal, @meg?

[Jaroslav Hron]

@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

[Jan Blechta]

@jxh's example works with W=FunctionSpace(mesh,'CG',1) but fails with W=FunctionSpace(Mesh(mesh),'CG',1). The first one collects the DOFs directly from vector when assembling rhs linear form while the latter one calls restrict_as_ufc_function(w, element, dolfin_cell, vertex_coordinates, ufc_cell); which must hence be broken also.

[Jan Blechta]

"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.

[Jan Blechta]

Incorrect code is not generated now. For correct reimplementation watch https://bitbucket.org/fenics-project/ffc/issue/69.

[Martin Sandve Alnæs]

Removing milestone: 1.6 (automated comment)