[firedrake] Possible bug with perturbed extruded mesh?
Tuomas Karna
tuomas.karna at gmail.com
Wed Apr 15 02:11:34 BST 2015
OK, had to pip install petsc and petsc4py to fix that.
Yes, the test works on deformed tetrahedral mesh:
Konsole output a == L: 0.99999977445 1.00000019589
a == L2: 0.999999999998 1.0
- Tuomas
On 04/14/2015 05:04 PM, Tuomas Karna wrote:
> 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