[firedrake] Possible bug with perturbed extruded mesh?
Tuomas Karna
tuomas.karna at gmail.com
Wed Apr 15 01:04:48 BST 2015
Hmm, just updated all the components, still seeing the same. I'm using
the mapdes firedrake branches for petsc and petsc4py. Could be linking
to old petsc4py or some other lib...
- Tuomas
On 04/14/2015 02:15 PM, Andrew McRae wrote:
> WTF?
>
> Is your PETSc (and everything else) up to date?
>
> I just opened an ipython session and typed
> from firedrake import *
> mesh = UnitCubeMesh(10, 10, 6)
>
> and all was fine.
>
> On 14 April 2015 at 21:51, Tuomas Karna <tuomas.karna at gmail.com
> <mailto:tuomas.karna at gmail.com>> wrote:
>
> Hi Andrew,
>
> I see, that could indeed be the reason.
>
> Tried UnitCubeMesh, results in an error:
>
> 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).
>>
>> 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 <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
More information about the firedrake
mailing list