The mesh is triangular, But the maximal element size (hmax) might differ from the value you Give to Gmsh. This is the h you should use in convergence studies.
Like that ?
mesh = dolfinx.UnitSquareMesh(MPI.COMM_WORLD, 10, 5)
tdim = mesh.topology.dim
num_cells = mesh.topology.index_map(tdim).size_local
h = dolfinx.cpp.mesh.h(mesh, tdim, range(num_cells))
Es = np.zeros(len(Ns_func), dtype=ScalarType)
list_rates_func = []
list_Es_func = []
dims_func = [1,2]
Ns = []
for n in range(len(dims_func)):
dim = dims_func[n]
for i in Ns:
Es[i] = errOnU(n_mesh=N, mesh_order=2, dim_func=dim)
list_Es_func.append(list(Es))
rates = np.log(Es[1:]/Es[:-1])/np.log(h[1:]/h[:-1])
list_rates_func.append(list(rates))
I think this is what @dokken was thinking about but the h
should of course be inside the loop which should be all about N
not dims_func
, if it’s convergence rates you’re looking for.
I would be surprised if this was the issue though. Since everything seems to work with P1 elements, I would more readily blame the variational formulation than the mesh size.
For the L2 norm I find : 1.92 for P1 and 3.45 for P2
while the theoretical results are : 2.00 for P1 and 3.00 for P2
For the H1 norm I find : 1.46 for P1 and 2.22 for P2
while the theoretical results are : 1.00 for P1 and 2.00 for P2
so i think nothing is working.
In fact I should be able to find the convergence rate and also plot the graphs as I showed you. I had chosen mesh_size = 2 instead of mesh_size = 1 because we observe better results
mesh = dolfinx.UnitSquareMesh(MPI.COMM_WORLD, 10, 5)
tdim = mesh.topology.dim
num_cells = mesh.topology.index_map(tdim).size_local
Es = np.zeros(len(Ns_func), dtype=ScalarType)
list_rates_func = []
list_Es_func = []
dims_func = [1,2]
Ns = []
for n in range(len(dims_func)):
dim = dims_func[n]
for i in Ns:
Es[i] = errOnU(n_mesh=N, mesh_order=2, dim_func=dim)
h = dolfinx.cpp.mesh.h(mesh, tdim, range(num_cells))
list_Es_func.append(list(Es))
rates = np.log(Es[1:]/Es[:-1])/np.log(h[1:]/h[:-1])
list_rates_func.append(list(rates))
If it’s not dim_func, how do you know it’s a polynomial of degree 1 or degree 2 that we want to see?
You should take some steps backwards.
You Need to first look at your PDE and see if there are theoretical convergence results for it in litterature.
Secondly. Start by looking at a problem were you can use the built in refinement routines of DOLFINx (first order, linear meshes).
Thirdly, simplify your problem, make sure that it works for a simple problem before progressing to your actual example.
Hi,
See B. Szabo, I. Babushka, Finite Element Analysis, 1991, pp. 69-70. They explain how to approximate the error in the energy norm without knowing the exact solution. For Poisson’s equation, the energy norm is exactly the H^1_0 norm.