### DG Quad elements do not converge

435
views
1
10 months ago by
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:

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],
# 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
- 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
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