← Back to team overview

dolfin team mailing list archive

Re: facet area in 1D

 

Looks like you need to check that you are not running in parallel in
that test since DOLFIN will otherwise try to partition the mesh which
doesn't work if you only have one cell.

--
Anders


On Fri, Nov 25, 2011 at 10:58:17AM +0100, Kristian Ølgaard wrote:
> On 24 November 2011 20:33, Anders Logg <logg@xxxxxxxxx> wrote:
> > 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.)
>
> I put the above script in
> test/unit/function/python/SpecialFunctions.py (not pushed to repo)
>
> but doing
>
> mpirun -np 3  python ./SpecialFunctions.py
>
> results in the error:
>
> Process 0: Number of global vertices: 2
> Process 0: Number of global cells: 1
> EE
> ======================================================================
> ERROR: testFacetArea (__main__.SpecialFunctions)
> ----------------------------------------------------------------------
> Traceback (most recent call last):
>   File "./SpecialFunctions.py", line 29, in testFacetArea
>     references = [(UnitInterval(1), 2, 2),]#\
>   File "/home/oelgaard/software/fenics/lib/python2.7/site-packages/dolfin/cpp.py",
> line 12959, in __init__
>     _cpp.UnitInterval_swiginit(self,_cpp.new_UnitInterval(*args))
> RuntimeError:
>
> *** -------------------------------------------------------------------------
> *** DOLFIN encountered an error. If you are not able to resolve this issue
> *** using the information listed below, you can ask for help at
> ***
> ***     https://answers.launchpad.net/dolfin
> ***
> *** Remember to include the error message listed below and, if possible,
> *** include a *minimal* running example to reproduce the error.
> ***
> *** -------------------------------------------------------------------------
> *** Error:  Unable to compute mesh partitioning using ParMETIS.
> *** Reason: ParMETIS cannot be used if a process has no cells (graph
> nodes). Use SCOTCH to perform partitioning instead.
> *** Where:  This error was encountered inside ParMETIS.cpp.
> *** -------------------------------------------------------------------------
>
>
> ----------------------------------------------------------------------
> Ran 1 test in 0.003s
>
> FAILED (errors=1)
>
> ======================================================================
> ERROR: testFacetArea (__main__.SpecialFunctions)
> ----------------------------------------------------------------------
> Traceback (most recent call last):
>   File "./SpecialFunctions.py", line 29, in testFacetArea
>     references = [(UnitInterval(1), 2, 2),]#\
>   File "/home/oelgaard/software/fenics/lib/python2.7/site-packages/dolfin/cpp.py",
> line 12959, in __init__
>     _cpp.UnitInterval_swiginit(self,_cpp.new_UnitInterval(*args))
> RuntimeError:
>
> *** -------------------------------------------------------------------------
> *** DOLFIN encountered an error. If you are not able to resolve this issue
> *** using the information listed below, you can ask for help at
> ***
> ***     https://answers.launchpad.net/dolfin
> ***
> *** Remember to include the error message listed below and, if possible,
> *** include a *minimal* running example to reproduce the error.
> ***
> *** -------------------------------------------------------------------------
> *** Error:  Unable to compute mesh partitioning using ParMETIS.
> *** Reason: ParMETIS cannot be used if a process has no cells (graph
> nodes). Use SCOTCH to perform partitioning instead.
> *** Where:  This error was encountered inside ParMETIS.cpp.
> *** -------------------------------------------------------------------------
>
>
> ----------------------------------------------------------------------
> Ran 1 test in 0.003s
>
> FAILED (errors=1)
>
>
> after which the process hangs.
>
> Kristian
>
> >


Follow ups

References