The below script shows some strange behaviour with 'auto' quadrature
degree and DG elements.
The script integrates du/dx over a unit square with u = the projection
of 'pi x' into some space.
For DG spaces of degree 1 and 2 and several choices of quadrature
orders, the result is either completely wrong or NaN.
Am I missing something here or should this be working better?
Also note that for CG elements the result of 'auto' quadrature degree is worse
than 0.
from dolfin import *
# Projection which takes form compiler parameters
def myproject(expr, space, params):
if 1: # Set to 0 to use the regular project and ignore params
U = Function(space)
u = TrialFunction(space)
v = TestFunction(space)
a = inner(u,v)*dx
L = inner(expr,v)*dx
problem = VariationalProblem(a, L, form_compiler_parameters=params)
problem.solve(U)
else:
U = project(expr, space)
return U
def test(fam, deg, qd):
ffc_opt = { 'representation': 'quadrature', 'quadrature_degree': qd }
n = 90
mesh = UnitSquare(n,n)
V = FunctionSpace(mesh, fam, deg)
u = myproject(3.14*triangle.x[0], V, ffc_opt)
a = u.dx(0)*dx
dudx = assemble(u.dx(0)*dx, form_compiler_parameters=ffc_opt)
print fam, deg, qd, dudx, norm(u.vector())
# Summary of the output is embedded in comments below:
print "Correct value is", 3.14
test("CG", 1, 0) # 3.14
test("CG", 1, 1) # 3.14
test("CG", 1, 2) # 3.14
test("CG", 1, 'auto') # 3.12 (worse than qd=0!)
test("DG", 2, 0) # NaN
test("DG", 2, 1) # NaN
test("DG", 2, 2) # NaN
test("DG", 2, 'auto') # 2.61 == 3.14/1.2
test("DG", 1, 0) # NaN
test("DG", 1, 1) # NaN
test("DG", 1, 2) # 3.14
test("DG", 1, 'auto') # almost zero
--
You received this bug notification because you are a member of DOLFIN
Team, which is subscribed to DOLFIN.
https://bugs.launchpad.net/bugs/785874
Title:
Projection of x is not accurate
Status in DOLFIN:
New
Bug description:
I've tested that projecting x works without the scaling bug that was
just fixed, using dimensions 1,2,3 and both DG and CG from 0 to 3
degrees. I print the max and min values of the vector of the
projection function, and the values are _close_ to 0 and 1 but not to
machine precision. The script is below.
There's up to 2.7% error in the 3D case. Is the projection form
integrated accurately enough? All but the DG0 function space should be
capable of representing x exactly. Not sure if this is a dolfin or ffc
bug.
from dolfin import *
def mcx(dim):
if dim == 1:
mesh = UnitInterval(20)
cell = interval
x = cell.x
if dim == 2:
mesh = UnitSquare(20, 20)
cell = triangle
x = cell.x[0]
if dim == 3:
mesh = UnitCube(20, 20, 20)
cell = tetrahedron
x = cell.x[0]
return mesh, cell, x
for dim in range(1, 4):
mesh, cell, x = mcx(dim)
minval, maxval = 1.0, 0.0
#print dim, "DG"
for degree in range(3):
V = FunctionSpace(mesh, "DG", degree)
u = project(x, V)
#print dim, degree, u.vector().min(), u.vector().max()
minval = min(minval, u.vector().min())
maxval = max(maxval, u.vector().max())
#print dim, "CG"
for degree in range(1, 4):
V = FunctionSpace(mesh, "CG", degree)
u = project(x, V)
#print dim, degree, u.vector().min(), u.vector().max()
minval = min(minval, u.vector().min())
maxval = max(maxval, u.vector().max())
print minval, maxval
_______________________________________________
Mailing list: https://launchpad.net/~dolfin
Post to : dolfin@lists.launchpad.net
Unsubscribe : https://launchpad.net/~dolfin
More help : https://help.launchpad.net/ListHelp
