Increasing the function-space order for Lagrange elements

variational formulation
1

Dear all,

I have a mixed function-space problem with complex numbers. Eventhough I increase the function-space order (p_h) the error of the numerical solution compared to the exact solution does not decrease at all for p_h > 3. It is somehow always stuck at p_h = 3 which should not (from other software and solutions for the same problem I know that). Here is how I define my function-space in my own problem:

Er = FiniteElement('P', triangle, p_h)
Ei = FiniteElement('P', triangle, p_h)
V = FunctionSpace(mesh,Ec)

Also, I define the exact solution function-space order as p_ex = p_h + 1. Both exact and numerical solutions are solved on the same mesh and the error is calculated by errornorm() function.

I was not able to write a simple code to show my problem here. But I found a topic closest to my problem. I have tried my best to do what they were describing but I was not able to reach to a better result in that problem either. Here is the topic that I think describes the same problem that I go through:

https://fenicsproject.org/qa/3795/solution-showing-linear-functions-higher-lagrange-elements/

The suggested solution is:
“The workaround I have been going through is to interpolate the higher order solution onto a finer linear space. If you need to write the solution on a finer mesh anyway, you could plot and write the finer mesh that’s an interpolation of the higher order solution.”

I don’t understand it. Do I need to create a new mesh to interpolate the high-order solution on-to? This is very confusing.

Here is the code that they have given and I have updated a little bit. I am trying to get a better result as their accepted solution suggests:

from dolfin import *
%matplotlib inline

mesh = UnitSquareMesh(3,3)

p = 3
q = 5

E = FiniteElement('CG', triangle, p)
V = FunctionSpace(mesh,E)

def boundary(x):
    return x[0] < DOLFIN_EPS or x[0] > 1.0 - DOLFIN_EPS

u0 = Constant(0.0)
bc = DirichletBC(V, u0, boundary)
u = TrialFunction(V)
v = TestFunction(V)
f = Expression("10*exp(-(pow(x[0] - 0.5, 4) + pow(x[1] - 0.5, 4)) / 0.02)", degree = q)
g = Expression("sin(5*x[0])", degree = q)
a = inner(grad(u), grad(v))*dx
L = f*v*dx + g*v*ds
u = Function(V)

solve(a == L, u, bc)

plot(u)

# From here on, it is my solution that looks like it will work
# but I want to project this solution (if it is actually the higher order solution) 
# to a finer mesh if I have to. I am not even sure from this point after

import numpy as np

W = FunctionSpace(mesh,'Lagrange', p+1)

w = interpolate(u, W)
print(w.vector()[:])
print(w.compute_vertex_values())

# Should I create a new mesh according to this higher-order solution? 
#I think I am doing what they suggest below, but I could not figure out 
# how to plot my new high-order solution to this new mesh also

gdim = mesh.geometry().dim()
dofmap = V.dofmap()
dofs = dofmap.dofs()
dofs_x = W.tabulate_dof_coordinates().reshape((-1, gdim))

# plot(?)

I hope I made it clear enough. I would like to plot the high-order solution function on my geometry in the end.

2

The code you are linking to is simply a visualization problem, you want to visualize higher order functions (matplotlib and Paraview interpolates the result into CG-1), which makes it tricky to visualize.

If you want more help on the topic, an example were the error stops reducing for higher p_h is needed.

3

Actually, I am using this mesh that you have coded for the most part:

But for FiniteElement(‘P’, triangle, p_h) for p_h > 3 it does not give me the results that I expect. I am assuming that my variational formulation is right. Even if it is wrong, that does not make sense to me that changing the order is not affecting the results at all. So I thought may be this is a common problem people are having with Lagrange elements when I found that topic.

4

I cannot really help you any further without an example. As your problem is defined on a cylinderical mesh, it might be that the linear approximation of the mesh geometry reduced convergence.

A suggestion for you is to try to make a manufactured solution for your PDE on a simple domain, see: Solving PDEs in Python - <br> The FEniCS Tutorial Volume I
for an example for Poisson problems, or Manufactured Solution for Navier Stokes - #2 by dokken for the Navier-Stokes problem with manufactured solution.

5

Dear dokken,

I was not able to describe my problem good enough but now I think I know what the problem is. The code below describes it a lot better this time.

from dolfin import *
from mshr import *
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline

# Define geometry
domain = Circle(Point(0, 0), 1, 30)

# Mesh the domain
mesh = generate_mesh(domain,10)

# Exact Solution Defined as a Function
def Exact_Solution(x,y):
    
    u_ex = pow(x,3)+2*y
    
    return u_ex

### LOWER ORDER EXACT SOLUTION ON MESH
p_ex = 1

F_1 = FunctionSpace(mesh, 'P', p_ex)

n_1 = F_1.dim()                                                                      
d_1 = mesh.geometry().dim()

e_1 = Function(F_1)

F_1_dof_coordinates = F_1.tabulate_dof_coordinates()
F_1_dof_coordinates.resize((n_1, d_1))
F_1_dof_x = F_1_dof_coordinates[:, 0]
F_1_dof_y = F_1_dof_coordinates[:, 1]

u_1_temp = Exact_Solution(F_1_dof_x,F_1_dof_y)
e_1.vector().set_local(u_1_temp)

# HIGHER ORDER EXACT SOLUTION ON MESH

F_2 = FunctionSpace(mesh, 'P', p_ex+1)

n_2 = F_2.dim()                                                                      
d_2 = mesh.geometry().dim()

e_2 = Function(F_2)

F_2_dof_coordinates = F_1.tabulate_dof_coordinates()
F_2_dof_coordinates.resize((n_1, d_1))
F_2_dof_x = F_2_dof_coordinates[:, 0]
F_2_dof_y = F_2_dof_coordinates[:, 1]

u_2_temp = Exact_Solution(F_2_dof_x,F_2_dof_y)
e_2.vector().set_local(u_2_temp)

plot(e_2-e_1)

# ERROR FOR CHANGING THE EXACT SOLUTION ORDER
print( 'Err e1 vs e2 =', errornorm(e_1, e_2, 'L2')/sqrt(assemble(pow(abs(e_2),2)*dx)))

# NUMERICAL MODEL

def Numeric_Solution(mesh,p_h):

    q = p_h
    V = FunctionSpace(mesh, 'P', p_h)
    u = interpolate(Expression('pow(x[0],3) + 2*x[1]', degree=q), V)
    
    return u

n_1 = Numeric_Solution(mesh, p_h=3)

# The error is bound to the order of the exact solution. Better the exact solution resolution, better the error analysis?
print( 'Err n1 vs e1 =', errornorm(e_1, n_1, 'L2')/sqrt(assemble(pow(abs(e_1),2)*dx)))
print( 'Err n1 vs e2 =', errornorm(e_2, n_1, 'L2')/sqrt(assemble(pow(abs(e_2),2)*dx)))

I think this is why my error is high. It is not realted to the Lagrange elements but to the method that I have to use to make an error analysis. My exact solution is a very complicated one. So I have to define it with a function. Then I solve it on the same mesh that I solve my numerical solution on. But when I compare the two solutions, I expect a lot lower error. Eventhough I increase the order of the exact solution, I never reach to the error that I can find in literature for my problem. Looks like I am limited with the exact solution order and mesh. I tried increasing the order of the exact solution and make a comparison yet it did not help. Still, the function-space order p_h > 3 does not change for p_ex = 8.

How can I project a numerical solution in a coarse mesh defined as a function in here to a finer mesh? If I am not doing something very fundamentally wrong ofcourse.