On Thu, Nov 24, 2011 at 01:43:20PM +0100, Kristian Ølgaard wrote: > 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
ok. Can you make it into a unit test? (Maybe you already did, I haven't checked.) -- Anders _______________________________________________ Mailing list: https://launchpad.net/~dolfin Post to : dolfin@lists.launchpad.net Unsubscribe : https://launchpad.net/~dolfin More help : https://help.launchpad.net/ListHelp
