[firedrake] Possible bug with perturbed extruded mesh?

[Andrew McRae]

Hi Tuomas,

It's possible that you're hitting a non-affine issue here.  At the moment,
the Jacobian, dx_i/dX_j, of the coordinate transformation from the
reference cell to the physical cell, x(X), is treated as being constant
over the cell.  For non-simplex cells (i.e. quadrilaterals, triangular
prisms, tetrahedra), this is only an approximation.  The Jacobian matrix is
used when replacing the measure (dx == |det J|*dX) and in derivatives
(d/dx_i == sum_j dX_j/dx_i d/dX_j).

I have some work in progress, hopefully landing in the next few weeks,
after which the Jacobian will be calculated separately at each quadrature
point.  This is currently spread across branches in different components of
Firedrake.  I'll try your code tomorrow to see if the results come out the
same (or at least much closer!).

Just checking, I assume that doing the problem on tetrahedra gives matching
results, even with a deformed mesh?  (If so, this suggests further that
it's a non-affine approximation issue).  I.e., replace the first four lines
of code with

from firedrake import *
mesh = UnitCubeMesh(10, 10, 6)  # mesh of tetrahedra, each of the 10*10*6
cuboids is split into 6 tets

and send it through the same deformation as before.  Replace 'ds_v + ds_t +
ds_b' by just ds in the declaration of L2.

Best,
Andrew

On 14 April 2015 at 19:29, Tuomas Karna <tuomas.karna at gmail.com> wrote:

> Hi all,
>
> Encountered this issue with extruded mesh where the coordinates have
> been deformed. The first form a==L just evaluates div(u), the second
> a==L2 is the same but div(u) integrated in parts. Thus the results
> should be equivalent, but that's not what I see:
>
> a == L: 0.999999886488 1.0000001092
> a == L2: -0.235329928101 3.47882106918
>
> Without the mesh deformation the two forms give the correct solution
> (1.0). Am I missing something?
>
> Cheers,
>
> Tuomas
>
> ---
>
> from firedrake import *
> mesh2d = UnitSquareMesh(10,10)
> layers = 6
> mesh = ExtrudedMesh(mesh2d, layers, -1.0/layers)
>
> P1 = FunctionSpace(mesh, 'CG', 1, vfamily='CG', vdegree=1)
> P2 = FunctionSpace(mesh, 'CG', 2, vfamily='CG', vdegree=1)
> P1v = VectorFunctionSpace(mesh, 'CG', 1, vfamily='CG', vdegree=1)
>
> # deform mesh
> scalar = Function(P1).interpolate(Expression('1.0+x[0]'))
> coords = mesh.coordinates
> coords.dat.data[:,2] *= scalar.dat.data[:]
>
> u = Function(P1v)
> w = Function(P2)
> u.interpolate(Expression(('x[0]', '0.0', '0.0')))
>
> tri = TrialFunction(P2)
> test = TestFunction(P2)
> normal = FacetNormal(mesh)
>
> a = test*tri*dx
> L = div(u)*test*dx
> L2 = -dot(u, grad(test))*dx + dot(u, normal)*test*(ds_v+ds_t+ds_b)
>
> solve(a == L, w)
> print 'a == L:', w.dat.data.min(), w.dat.data.max()
>
> solve(a == L2, w)
> print 'a == L2:', w.dat.data.min(), w.dat.data.max()
>
>
> _______________________________________________
> firedrake mailing list
> firedrake at imperial.ac.uk
> https://mailman.ic.ac.uk/mailman/listinfo/firedrake
>
-------------- next part --------------
HTML attachment scrubbed and removed