[firedrake] Possible bug with perturbed extruded mesh?

[Tuomas Karna]

Hi Andrew,

I see, that could indeed be the reason.

Tried UnitCubeMesh, results in an error:

Konsole output File "test_integration.py", line 77, in <module>
     mesh = UnitCubeMesh(10, 10, 6)
   File 
"/home/tuomas/.local/lib/python2.7/site-packages/firedrake/utility_meshes.py", 
line 508, in UnitCubeMesh
     return CubeMesh(nx, ny, nz, 1, reorder=reorder)
   File 
"/home/tuomas/.local/lib/python2.7/site-packages/firedrake/utility_meshes.py", 
line 488, in CubeMesh
     return BoxMesh(nx, ny, nz, L, L, L, reorder=reorder)
   File 
"/home/tuomas/.local/lib/python2.7/site-packages/firedrake/utility_meshes.py", 
line 440, in BoxMesh
     plex = PETSc.DMPlex().generate(boundary)
   File "PETSc/DMPlex.pyx", line 438, in petsc4py.PETSc.DMPlex.generate 
(src/petsc4py.PETSc.c:215426)
petsc4py.PETSc.Error: error code 77
[0] DMPlexGenerate() line 1080 in 
/home/tuomas/sources/petsc/firedrake/src/dm/impls/plex/plexgenerate.c
[0] DMPlexGenerate_CTetgen() line 834 in 
/home/tuomas/sources/petsc/firedrake/src/dm/impls/plex/plexgenerate.c
[0] TetGenTetrahedralize() line 21483 in 
/home/tuomas/sources/petsc/firedrake/arch-linux2-c-debug/externalpackages/ctetgen/ctetgen.c
[0] TetGenMeshDelaunizeVertices() line 12113 in 
/home/tuomas/sources/petsc/firedrake/arch-linux2-c-debug/externalpackages/ctetgen/ctetgen.c
[0] TetGenMeshDelaunayIncrFlip() line 12046 in 
/home/tuomas/sources/petsc/firedrake/arch-linux2-c-debug/externalpackages/ctetgen/ctetgen.c
[0] TetGenMeshInsertVertexBW() line 11321 in 
/home/tuomas/sources/petsc/firedrake/arch-linux2-c-debug/externalpackages/ctetgen/ctetgen.c
[0] TetGenMeshPreciseLocate() line 5781 in 
/home/tuomas/sources/petsc/firedrake/arch-linux2-c-debug/externalpackages/ctetgen/ctetgen.c
[0] Petsc has generated inconsistent data
[0] This is wrong


On 04/14/2015 12:08 PM, Andrew McRae wrote:
> 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).
>
> Konsole output 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 
> <mailto: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 <mailto:firedrake at imperial.ac.uk>
>     https://mailman.ic.ac.uk/mailman/listinfo/firedrake
>
>
>
>
> _______________________________________________
> firedrake mailing list
> firedrake at imperial.ac.uk
> https://mailman.ic.ac.uk/mailman/listinfo/firedrake

-------------- next part --------------
HTML attachment scrubbed and removed