On 23 November 2011 10:54, Martin Sandve Alnæs <marti...@simula.no> wrote:
> On 23 November 2011 09:37, Kristian Ølgaard <k.b.oelga...@gmail.com> wrote:
>> On 22 November 2011 21:15, Kristian Ølgaard <k.b.oelga...@gmail.com> wrote:
>>> On 22 November 2011 21:07, Anders Logg <l...@simula.no> wrote:
>>>> On Tue, Nov 22, 2011 at 07:34:38PM +0000, Garth N. Wells wrote:
>>>>> On 22 November 2011 19:12, Kristian Ølgaard <k.b.oelga...@gmail.com>
>>>>> wrote:
>>>>> > Hi,
>>>>> >
>>>>> > The IntervalCell::facet_area in DOLFIN currently returns 0.0, while
>>>>> > the facet determinant (a scale factor) 'det = 1.0;'
>>>>> > in the generated code for the tabulate_tensor function in 1D (facet
>>>>> > integrals).
>>>>> >
>>>>> > Is 0.0 the correct return value, or should it return 1.0?
>>>>> >
>>>>>
>>>>> I would think 1.0.
>>>
>>> That's what I'm leaning towards too.
>>>
>>>> Facet is the thing of dimension one lower than the cell dimension so
>>>> in this case 1 - 1 so it would mean the area of an endpoint. So I
>>>> believe 0 is correct.
>>>
>>> I see your 'point' :), but is it consistent with how we evaluate
>>> integrals on endpoints?
>>
>> Running this script:
>>
>> from dolfin import *
>>
>> cells = [interval, triangle, tetrahedron]
>> meshes = [UnitInterval(1), UnitSquare(1,1), UnitCube(1,1,1)]
>>
>> for cell, mesh in zip(cells, meshes):
>
> btw, you can do
> cell = mesh.ufl_cell()
>
>> c = Constant(1, cell)
>> print
>> print cell
>> print "volume: ", assemble(c*dx, mesh=mesh)
>> print "surface: ", assemble(c*ds, mesh=mesh)
>>
>> produces:
>>
>> <interval cell in R1>
>> volume: 1.0
>> surface: 2.0
>>
>> <triangle cell in R2>
>> volume: 1.0
>> surface: 4.0
>>
>> <tetrahedron cell in R3>
>> volume: 1.0
>> surface: 6.0
>>
>> If we want surface = 0.0 for the interval, 'det' must be set to zero
>> in the generated code.
>> However, doing so will mean that f*ds etc. becomes zero too, which is
>> not very useful.
>>
>> Kristian
>
> As Kristian shows here, I think defining it to 1.0 seems more useful,
> and consistent with the definition of 1.0*ds.
I've pushed the change to DOLFIN, the output of the extended script:
from dolfin import *
for mesh in [UnitInterval(1), UnitSquare(1,1), UnitCube(1,1,1)]:
c = Constant(1, mesh.ufl_cell())
c0 = mesh.ufl_cell().facet_area
c1 = FacetArea(mesh)
print
print mesh.ufl_cell()
print "volume: ", assemble(c*dx, mesh=mesh)
print "surface: ", assemble(c*ds, mesh=mesh)
print "facet_area ufl/ffc: ", assemble(c0*ds, mesh=mesh)
print "FacetArea dolfin: ", assemble(c1*ds, mesh=mesh)
is
<interval cell in R1>
volume: 1.0
surface: 2.0
facet_area ufl/ffc: 2.0
FacetArea dolfin: 2.0
<triangle cell in R2>
volume: 1.0
surface: 4.0
facet_area ufl/ffc: 4.0
FacetArea dolfin: 4.0
<tetrahedron cell in R3>
volume: 1.0
surface: 6.0
facet_area ufl/ffc: 3.0
FacetArea dolfin: 3.0
Kristian
> Martin
>
_______________________________________________
Mailing list: https://launchpad.net/~dolfin
Post to : dolfin@lists.launchpad.net
Unsubscribe : https://launchpad.net/~dolfin
More help : https://help.launchpad.net/ListHelp
