### DG Quad elements do not converge

435

views

1

I'm trying to solve DG problems using the new DGQuad elements, however I can't seem to get them to converge to a known solution.

I've constructed a basic manufactured solution test for the linear Poisson equation,

editing the demo code slightly.

I also assemble the cell volumes into a DG0 function since CellSize and CellVolume are currently not working with DGQuad elements.

In the following code, turn off/on the UnitSquareMesh and UnitQuadMesh to demonstrate the effect.

If anyone can spot something stupid that I've done, that would be wonderful.

Edited code for FEniCS master 3-Apr-2018:

I've constructed a basic manufactured solution test for the linear Poisson equation,

editing the demo code slightly.

I also assemble the cell volumes into a DG0 function since CellSize and CellVolume are currently not working with DGQuad elements.

In the following code, turn off/on the UnitSquareMesh and UnitQuadMesh to demonstrate the effect.

If anyone can spot something stupid that I've done, that would be wonderful.

Edited code for FEniCS master 3-Apr-2018:

```
from dolfin import *
import numpy as np
parameters["ghost_mode"] = "shared_facet"
run_count = 0
ele_ns = [4, 8, 16, 32, 64]
errorl2 = np.zeros(len(ele_ns))
errorh1 = np.zeros(len(ele_ns))
hsizes = np.zeros(len(ele_ns))
for ele_n in ele_ns:
# Create mesh and define function space
mesh = RectangleMesh.create(MPI.comm_world,
[Point(0, 0), Point(1, 1)], [ele_n, ele_n],
CellType.Type.quadrilateral)
# mesh = UnitSquareMesh(ele_n, ele_n)
V = FunctionSpace(mesh, 'DG', 2)
# Define test and trial functions
u = TrialFunction(V)
v = TestFunction(V)
# Assemble cell volumes into a DG0 function since quad elements do not support
# CellSize() nor CellVolume()
DG0 = FunctionSpace(mesh, "DG", 0)
cell_volume = Function(DG0)
v0 = TestFunction(DG0)
cell_volume.vector()[:] = assemble(v0*dx)
# # Define normal vector and mesh size
n = FacetNormal(mesh)
h = cell_volume/FacetArea(mesh)
h_avg = (h('+') + h('-'))/2
# Define the source term f, Dirichlet term u0 and Neumann term g
f = Expression('2*pi*pi*sin(pi*x[0])*sin(pi*x[1])', element=V.ufl_element())
u0 = Expression("sin(pi*x[0])*sin(pi*x[1])", element=V.ufl_element())
# Mark facets of the mesh
boundaries = MeshFunction('size_t', mesh, 1, 0)
AutoSubDomain(lambda x, on: on).mark(boundaries, 1)
# Define outer surface measure aware of Dirichlet and Neumann boundaries
ds = Measure('ds', domain=mesh, subdomain_data=boundaries)
# Define parameters
alpha = 4.0
gamma = 8.0
# Define variational problem
a = dot(grad(v), grad(u))*dx \
- dot(avg(grad(v)), jump(u, n))*dS \
- dot(jump(v, n), avg(grad(u)))*dS \
+ alpha/h_avg*dot(jump(v, n), jump(u, n))*dS \
- dot(grad(v), u*n)*ds(1) \
- dot(v*n, grad(u))*ds(1) \
+ (gamma/h)*v*u*ds(1)
L = v*f*dx - u0*dot(grad(v), n)*ds(1) + (gamma/h)*u0*v*ds(1)
# Compute solution
u = Function(V)
solve(a == L, u)
# Assemble l2 error
errorl2[run_count] = assemble((u0 - u)**2*dx)**0.5
hsizes[run_count] = mesh.hmax()
run_count += 1
# Compute convergence rates. L2 SIPG convergence rate should be ~ p+1
if MPI.rank(mesh.mpi_comm()) == 0:
print(np.log(errorl2[0:-1]/errorl2[1:])/np.log(hsizes[0:-1]/hsizes[1:]))
```

Community: FEniCS Project

1

I have seen it, looks critical. https://bitbucket.org/fenics-project/dolfin/pull-requests/406/deprecate-cellsize-and-introduce/diff#comment-45866214

written
10 months ago by
Jan Blechta

1

This is currently being discussed at https://bitbucket.org/fenics-project/ffc/pull-requests/91

written
9 months ago by
Jan Blechta

1

Does it work with current master branches?

written
4 months ago by
Jan Blechta

1

I edited the code to work with the master branch(es). Seems to work correctly and I see optimal convergence rates on quads.

written
4 months ago by
Nate

Please login to add an answer/comment or follow this question.