Projection from fine low order basis function to coarse high order basis function

Hi, Artur Palha

    I think FEniCS has done what you expected. The plots always use lines to show the curves.

If you try to check the order of errors between the exact function and the projected function, the results agree with the analysis, i.e., the projection happened and provided
correct values.
Please see the code below

#----------------------------

from dolfin import *

# define the expression of the function to interpolate
F_analytical = Expression("sin(2*pi*x[0])")

# define the coarse and fine meshes
nx_list = [20,40,80]

for nx in nx_list:

    mesh_fine = UnitIntervalMesh(nx)
    mesh_coarse = UnitIntervalMesh(nx/2)

    # define the low and high order function spaces
    V_low_order = FunctionSpace(mesh_fine,"Lagrange",1)
    V_high_order = FunctionSpace(mesh_coarse,"Lagrange",2)

    # project the function into the low order basis function space
    F_low_order = project(F_analytical,V_low_order)

    # project the function into the high order basis function space
    F_high_order = project(F_analytical,V_high_order)

    l2error = errornorm(F_analytical, F_high_order, 'L2', degree_rise=3)

    print l2error

#----------------------------

On Wed, 10 Apr 2013 09:11:16 -0000
Artur Palha <email address hidden> wrote:
> New question #226361 on DOLFIN:
> https://answers.launchpad.net/dolfin/+question/226361
>
> I have data in a fine mesh which I triangulate and interpolate with
> linear basis functions. Then I want to represent this data with a
> coarser mesh but with higher order basis functions. Typically the
> number of degrees of freedom in my coarser higher order mesh is equal
> or smaller than the one in the finer low order mesh.
>
> If I do a Least Squares finite element approach to this problem I
> should get something like this:
>
> <u_l,v> = <u_h,v>
>
> Where u_l is the finite element representation of my function in the
> fine low order basis and u_h is the representation of my function in
> the coarse high order basis (the least squares approximation I wish
> to find) and v are the high order basis functions.
>
> But this is what, I think, the function projection does. So, I tried
> it and it did not work. I get a piecewise linear representation of
> the fine mesh but in a coarser mesh. So no use for the high order
> basis.
>
> This is an example in 1D of what I am talking about. I was expecting
> F_high_from_low_order to be high order, but it is not. Can anyone
> explain why this happens?
>
> Thank you
>
> from dolfin import *
>
> # define the expression of the function to interpolate
> F_analytical = Expression("sin(2*pi*x[0])")

You might want to increase Expression degree form default 1.
F_analytical = Expression("sin(2*pi*x[0])", degree=2)
It's used by UFL to estimate necessary quadrature rule. UFL can't see
compiled expression C++ code to analyse it automatically.

>
> # define the coarse and fine meshes
> mesh_fine = UnitIntervalMesh(12)
> mesh_coarse = UnitIntervalMesh(4)
>
> # define the low and high order function spaces
> V_low_order = FunctionSpace(mesh_fine,"Lagrange",1)
> V_high_order = FunctionSpace(mesh_coarse,"Lagrange",3)
>
> # project the function into the low order basis function space
> F_low_order = project(F_analytical,V_low_order)
>
> # project the function into the high order basis function space
> F_high_order = project(F_analytical,V_high_order)
>
> # project from high order to low order to plot with higher resolution
> F_high_order_plot = project(F_high_order,V_low_order)
>
> # project from low order to high order
> F_high_from_low_order = project(F_low_order,V_high_order)
>
> # project this new high order into lower order to plot with higher
> resolution F_high_from_low_order_plot =
> project(F_high_from_low_order,V_low_order)
>
> # plot everything
> plot(F_low_order)
> plot(F_high_order_plot)
> plot(F_high_from_low_order_plot)

Here's a primary problem: plots are always made with CG1 functions
interpolated to vertex values as well as PVD output.

Try this:

from dolfin import *
mesh = UnitIntervalMesh(1)
V = FunctionSpace(mesh, 'CG', 2)
u = Function(V)
u.vector()[0] = -100.
u.vector()[1] = 10.
u.vector()[2] = 100.
plot(u) # plots a linear function from 10 to 100 but
u(0.) # 10
u(1.) # 100
u(0.5) # -100

>
> interactive()
>
>

Thank you for you replies Jan and Huadong, but I think you did not understand my point.

I know that if I have a function represented with high order basis functions and plot it, only the values at the vertices are plotted, therefore the plot "seems" linear.

For that reason, whenever I want to plot a high resolution of a high order function I project it (or interpolate it) into a finer mesh. This is plotting issues and this is not the issue.

The issue is (and I took the code of Huadong and changed it to express my problem:

The code is exactly as yours. The only differences are that I plot the solutions and print the error of all solutions.

There are 3 solutions:

1- linear basis functions on a fine mesh (interpolated directly from the analytical function)
2- quadratic basis functions on a coarse mesh (inteprolated directly from the analytical function)
3- quadratic basis functions on a coarse mesh (projected from 1 (linears on fine mesh) into 2 (quadratics on coarse mesh))

Now, if I plot all of them just with plot I see that 1 has a finer detail but 2 and 3 not. This is natural because the higher order functions only plot the values at the vertices. For that reason when I plot higher order functions I project them to a finer lower order mesh (just for plotting).

Ok, now if you look at the plots generated by the code bellow you see that:

1- PLOT 1 has the low order interpolation of the analytical function. (this is correct)

2- PLOT 4 has the high order interpolation of the analytical function, plotted on a higher resolution mesh. So that we can see the higher detail. (this is correct)

3- PLOT 3 has the high order interpolation of the analytical function. Direct plot, which shows only the vertices. (this is correct)

4- PLOT 2 has the projection of the low order basis (fine mesh) into the high order basis (coarse mesh), plotted on a higher resolution mesh. This plot should be very similar to PLOT 4 (refered to in 2). But instead is like PLOT 3. This high order representation is linear, not quadratic.

Can you explain this?

This line:
F_high_order_from_low_order = project(F_low_order,V_high_order)

is just a Least Squares approximation of the low order represention into a higher order representation. Why does the result forget about the interior points?

Thank you once again.

from dolfin import *

# define the expression of the function to interpolate
F_analytical = Expression("sin(2*pi*x[0])",degree=2)

# define the coarse and fine meshes
nx_list = [10,20,40]

for nx in nx_list:

    mesh_fine = UnitIntervalMesh(nx)
    mesh_coarse = UnitIntervalMesh(nx/2)

    # define the low and high order function spaces
    V_low_order = FunctionSpace(mesh_fine,"Lagrange",1)
    V_high_order = FunctionSpace(mesh_coarse,"Lagrange",2)

    # project the function into the low order basis function space
    F_low_order = project(F_analytical,V_low_order)

    # project the function into the high order basis function space
    F_high_order = project(F_analytical,V_high_order)

    # project low order to high order
    F_high_order_from_low_order = project(F_low_order,V_high_order)

    l2error_high = errornorm(F_analytical, F_high_order, 'L2', degree_rise=4)
    l2error_low = errornorm(F_analytical, F_low_order, 'L2', degree_rise=4)
    l2error_high_from_low = errornorm(F_analytical, F_high_order_from_low_order, 'L2', degree_rise=3)

    plot(F_low_order,title='Low order basis')
    plot(project(F_high_order_from_low_order,V_low_order),title='High order (analytical) high resolution plot') # project the (exact) higher order to low order for plotting
    plot(F_high_order_from_low_order,title='High order (projected from low order) high resolution plot') # plot the high order projected from the low order

    print l2error_high, l2error_low, l2error_high_from_low

    interactive()

Hi Artur

    I have read your post, and I also run the above code. I still think everything is correct.

I believe you question is the following
"This high order representation is linear, not quadratic. "

Since I am in Department of mathematics, please allow me to do some analysis.

let F_low_order denotes the projection of exact F into fine mesh with low basis, and let h be the meshsize, so the L2-norm error of these two is O(h^2), that is we should have second order convergence.

Now, we project this F_low_order into a coarse mesh with high order basis, for example second order,
 let us use F_high_order_from_low_order denotes the projection of F_low_order, so the L2-norm error of these two is O((2h)^3),

Finally, we want to know how close F_high_order_from_low_order to the exact F, from the above analysis, it should be O(h^2+(2h)^3). Therefore, we can only expect a second order convergence.

Hope I can provide some real help. I am sorry if I misunderstand your question.

Huadong,

Thanks for your time.

I tried something similar to this with my own finite element code (in 2D on quads). And, for the case of these mesh relations I get higher order convergence. And I think the higher order convergence is correct.

Look at the simple example of a parabola: f(x) = -x*x + x.

With two linear elements on the unit interval, the F_low_order will be simply: f_0 = 0, f_1 = 1, f_2 =0. So, a linear hat. If we try to find the Least Squares fit of this hat with a parabola. The result is a parabola that passes very close to the nodes of the linear hat. Since those nodes are the exact nodes (or very close to it), the error is zero (or very close to zero).

Now, if you try that simple example with dolfin you will get: 0. That is f(x) = 0 for the quad. Which makes no sense.

The Least Square Finite Element approach (which is the same as the projection function) gives the weak formulation:

<f_low_order, v_high_order> = <f_high_from_low, v_high_order>

If you assemble this in dolfin for the simple mesh I give, the left hand side gives zero. Which for me makes no sense.

What is exactly being computed in <f_low_order, v_high_order>? Does anybody know exactly how this is done when the meshes are not the same?

Because what we have is:

   element 0 element1
|-----------------|-----------------| ----> fine mesh (low order)

              element 0
|-----------------------------------| ----> coarse mesh (high order)

Where are the quadrature points when we want to compute: <f_low_order, v_high_order>?

Moreover, if I try to assemble the following inner product matrix:

<v_test_low, v_trial_high>

By using:

L = inner(v_test_low,v_trial_high)*dx
assemble(L)

Then I get an error:

*** -------------------------------------------------------------------------
*** Error: Unable to extract mesh from form.
*** Reason: Non-matching meshes for function spaces.
*** Where: This error was encountered inside Form.cpp.
*** Process: 0
*** -------------------------------------------------------------------------

Which makes sense since the meshes are non-matching. (it might be a functionality not implemented in dolfin). Now, what does not makes sense is that this does not work but project works. Is project interpolating to a common mesh (the coarse one) and then projecting? That is why I get the low order on the higher order mesh?

Hi Artur

I try to compare the project of FEniCS and our own L2 prjection,
i.e., we can project the function F into the FEM space in L2-norm.

Please see the code below, and you can see that the command
project in FEniCS is exactly the L2 projection.

from dolfin import *

# define the expression of the function to interpolate
F_analytical = Expression("sin(2*pi*x[0])",degree=5)

# define the coarse and fine meshes
nx_list = [20,40,80]

for nx in nx_list:

    mesh_fine = UnitIntervalMesh(nx)
    mesh_coarse = UnitIntervalMesh(nx/2)

    # define the low and high order function spaces
    V_low_order = FunctionSpace(mesh_fine,"Lagrange",1)
    V_high_order = FunctionSpace(mesh_coarse,"Lagrange",2)

    # project the function into the low order basis function space (FEniCS's projection)
    F_low_order = project(F_analytical,V_low_order)

    # project the function into the low order basis function space (L2 projection)
    u = TrialFunction(V_low_order)
    v = TestFunction(V_low_order)
    a = u*v*dx
    L = F_analytical*v*dx
    u = Function(V_low_order)
    solve(a == L, u)

    l2error_l2projection = errornorm(F_analytical, u, 'L2', degree_rise=2)
    l2error_low = errornorm(F_analytical, F_low_order, 'L2', degree_rise=2)

    print l2error_l2projection, l2error_low

Huadong,

But what you show is just the projection of an analytical function into a FEM space. That works fine.

What I say is that the projection between FEM spaces defined in two distinct meshes is working in a strange manner for me.

This the the simplest code I can come up with. Why is F_high_from_low linear?

from dolfin import *

F_exact = Expression('-x[0]*x[0]+2*x[0] + 1.0')

mesh_low = UnitIntervalMesh(2)
mesh_high = UnitIntervalMesh(1)

V_low = FunctionSpace(mesh_low,'Lagrange',1)
V_high = FunctionSpace(mesh_high,'Lagrange',2)

v_trial_low = TrialFunction(V_low)
v_test_low = TestFunction(V_low)

v_trial_high = TrialFunction(V_high)
v_test_high = TestFunction(V_high)

F_low = project(F_exact,V_low)
F_high = project(F_exact,V_high)
F_high_from_low = project(F_low,V_high)

plot(F_low)
plot(interpolate(F_high,V_low))
plot(interpolate(F_high_from_low,V_low))

interactive()

Hi Artur

    I think I understand your question now.

    Yes, you are correct. Thank you for pointing this to me.

    I also tried another function f(x)=sin(2*x), similar results are obtained. I do not how the projection is done with a discrete function. There maybe something wrong in the projection between different FEM spaces, I have no idea.

    Let us hope someone know the details.

On Thu, 11 Apr 2013 14:51:05 -0000
Huadong GAO <email address hidden> wrote:
> Question #226361 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/226361
>
> Huadong GAO posted a new comment:
> Hi Artur
>
> I think I understand your question now.
>
> Yes, you are correct. Thank you for pointing this to me.
>
> I also tried another function f(x)=sin(2*x), similar results are
> obtained. I do not how the projection is done with a discrete
> function. There maybe something wrong in the projection between
> different FEM spaces, I have no idea.

Let u some function. Than projection of u to space V is u_proj
solving linear variational problem: inner(u_proj, v) == inner(u, v) for
all v from V (subject to bcs if specified).

There's not something wrong between different spaces. At first there
are not two spaces involved. You have only one space V and
it suffices that generic function u can be evaluated at needed
quadrature points. It does not even need to be from some function space.

Only exception is that projection, interpolation and generally function
evaluation does not work on non-matching meshes in parallel.

Jan

>
> Let us hope someone know the details.
>

Jan Blechta,

Thank you for your time.

I cannot understand what you want to say. Do you mean that the projection does not work if we use non matching meshes?

If so, why does it work the other way around? That is, when I project the high order basis on a coarse mesh, into a low order basis on a fine mesh?

Do you know how the inner products are computed for non-matching meshes?

Moreover, how can I compute the high order basis functions on a set of points? Because if this was possible, the problem could be solved.

Thank you

-artur palha

On Thu, 11 Apr 2013 15:35:58 -0000
Artur Palha <email address hidden> wrote:
> Question #226361 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/226361
>
> Status: Answered => Open
>
> Artur Palha is still having a problem:
> Jan Blechta,
>
> Thank you for your time.
>
> I cannot understand what you want to say. Do you mean that the
> projection does not work if we use non matching meshes?

No. I said that it does not work on non-matching meshes in PARALLEL. It
works serially.

>
> If so, why does it work the other way around? That is, when I project
> the high order basis on a coarse mesh, into a low order basis on a
> fine mesh?
>
> Do you know how the inner products are computed for non-matching
> meshes?

Sorry, I forgot *dx. These are normal UFL forms:
  lhs = inner(u_proj, v)*dx
  rhs = inner(u, v)*dx

lhs is bilinear form with
u_proj, v = TrialFunction(V), TestFunction(V)

rhs is linear form with v = TestFunction(V) and u being GenericFunction
which can be evaluated on needed quadrature points. Therefore
it suffices if u can be evaluated on V's mesh.

>
> Moreover, how can I compute the high order basis functions on a set of
> points? Because if this was possible, the problem could be solved.

Try something of:
BasisFunction class
FiniteElement.evaluate_basis()
FiniteElement.evaluate_basis_all()
FiniteElement.evaluate_basis_derivatives()
FiniteElement.evaluate_basis_derivatives_all()

Check documentation or try grep demo dir.

Jan

>
> Thank you
>
> -artur palha
>

This might be relevant:

    https://bugs.launchpad.net/dolfin/+bug/901069

Thanks to all you guys, especially Garth for the very handy link.

The problem they report is my problem (the only thing is that I use high order basis on the coarse mesh). I used the inneficient fix that was given there. I put it here. This solver the problem I have. The only thing is that it is inneficient to loop over all degrees of freedom in python. Isn't there a way to do this efficiently?

Once again thank you all.

mesh_coarse = UnitIntervalMesh(1)
V_coarse = FunctionSpace(mesh_coarse, "CG", 2)
mesh_fine = UnitIntervalMesh(2)
V_fine = FunctionSpace(mesh_fine, "CG", 1)

f = Expression('-x[0]*x[0]+2*x[0] + 1.0')
f_fine = interpolate(f, V_fine)
f_coarse = interpolate(f, V_coarse)

u = TrialFunction(V_coarse)
v = TestFunction(V_coarse)
lhs = u*v*dx
A = assemble(lhs)

b_correct = Vector(3)
for i in range(3):
    basis = Function(V_coarse)
    basis.vector()[i] = 1.0
    b_correct[i] = assemble(f_fine*basis*dx, mesh = mesh_fine)
f_correct = Function(V_coarse)
solve(A, f_correct.vector(), b_correct)

plot(interpolate(f_correct,V_fine))
interactive()