dolfin team mailing list archive

[Anders Logg <logg@xxxxxxxxx>]

On Thu, Nov 24, 2011 at 01:43:20PM +0100, Kristian Ølgaard wrote:
> On 23 November 2011 10:54, Martin Sandve Alnæs <martinal@xxxxxxxxx> wrote:
> > On 23 November 2011 09:37, Kristian Ølgaard <k.b.oelgaard@xxxxxxxxx> wrote:
> >> On 22 November 2011 21:15, Kristian Ølgaard <k.b.oelgaard@xxxxxxxxx> wrote:
> >>> On 22 November 2011 21:07, Anders Logg <logg@xxxxxxxxx> 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.oelgaard@xxxxxxxxx> 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

ok.

Can you make it into a unit test? (Maybe you already did, I haven't checked.)

--
Anders