dolfin team mailing list archive
-
dolfin team
-
Mailing list archive
-
Message #23577
Re: [Bug 785874] Re: Projection of x is not accurate
Ok, so NaN was because I manually specified too low degree.
So we're left with 'auto' which gives a big error:
0.00666666666667 0.993333333333 8.77744057992e-16
This is probably caused by ffc not considering
spatial coordinates in the expression. Fair enough.
So auto probably chooses 2 for the bilinear form
and 1 for the linear form, keeping the matrix nonsingular
but reducing the accuracy of the right hand side.
So the only bug here is that FFC could look at
spatial coordinates in the expression for better
autodetection, which would make project() of
expressions in PyDOLFIN more robust.
I'll close this bug and file a blueprint with this point.
Martin
On 30 May 2011 22:22, Garth Wells <785874@xxxxxxxxxxxxxxxxxx> wrote:
> The NaN is not caused by FFC autodetecting the degree.
>
> The numerical error in some cases is due to FFC not increasing the
> degree when the spatial coordinate is used as a coefficient function.
>
> --
> You received this bug notification because you are a member of DOLFIN
> Core 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
>
--
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
Follow ups
References