dolfin team mailing list archive

[Johan Hake <hake@xxxxxxxxx>]

On Wednesday 19 August 2009 23:26:56 Garth N. Wells wrote:
> Johan Hake wrote:
> > On Wednesday 19 August 2009 15:48:43 Garth N. Wells wrote:
> >> Johan Hake wrote:
> >>> On Wednesday 19 August 2009 15:12:21 Garth N. Wells wrote:
> >>>> Johan Hake wrote:
> >>>>> On Wednesday 19 August 2009 14:57:20 Garth N. Wells wrote:
> >>>>>> Johan Hake wrote:
> >>>>>>> On Wednesday 19 August 2009 14:21:12 Garth N. Wells wrote:
> >>>>>>>> Garth N. Wells wrote:
> >>>>>>>>> Kent Andre wrote:
> >>>>>>>>>>> Not sure what we have to do actually. I suppose I want get any
> >>>>>>>>>>> clue from the
> >>>>>>>>>>>   demo/function/restriction
> >>>>>>>>>>>
> >>>>>>>>>>> demo?
> >>>>>>>>>>> I see that Kent has modified the dolfin.FunctionSpace to
> >>>>>>>>>>> reflect the present (broken) implementation of restriction,
> >>>>>>>>>>> anything in that line? Should such an implementation has been
> >>>>>>>>>>> done in FunctionSpaceBase instead of FunctionSpace? Kent?
> >>>>>>>>>>>
> >>>>>>>>>>> Johan
> >>>>>>>>>>
> >>>>>>>>>> This restriction is different in the sense that the restriction
> >>>>>>>>>> is made by the standard element on a sub mesh.
> >>>>>>>>>> The restriction in focus now is on the complete mesh but for
> >>>>>>>>>> only a part of the element.
> >>>>>>>>>> Combined, these two features will be very powerful.  I guess the
> >>>>>>>>>> sub-mesh restriction stuff is broken due to the parallel mesh.
> >>>>>>>>>
> >>>>>>>>> I've attached an example solver. To make it work for P^k with k >
> >>>>>>>>> 2 we need to restrict some functions to cell facets. Using FFC,
> >>>>>>>>> we do
> >>>>>>>>>
> >>>>>>>>>     scalar  = FiniteElement("Lagrange", "triangle", 3)
> >>>>>>>>>     dscalar = FiniteElement("Discontinuous Lagrange", "triangle",
> >>>>>>>>> 3)
> >>>>>>>>>
> >>>>>>>>>     # This syntax restricts the space 'scalar' to cell facets.
> >>>>>>>>>     mixed = scalar["facet"] + dscalar
> >>>>>>>>>
> >>>>>>>>> In PyDOLFIN we could do
> >>>>>>>>>
> >>>>>>>>>     V_dg = FunctionSpace(mesh, "DG", 3)
> >>>>>>>>>     V_cg = FunctionSpace(mesh, "CG", 3, "facet")
> >>>>>>>>>
> >>>>>>>>>     mixed = V_cg + V_dg
> >>>>>>>>
> >>>>>>>> I see now that we need to have the above option since the
> >>>>>>>> everything (e.g. the dof map) is created and compiled inside
> >>>>>>>> functionspace.py.
> >>>>>>>
> >>>>>>> We could implement it in FunctionSpaceBase, as an operator,
> >>>>>>> returning a new FunctionSpaceBase, with the restricted element as
> >>>>>>> arguemnt.
> >>>>>>>
> >>>>>>> However my concerns go to a sub space of a restricted mixed space
> >>>>>>> (if such exists) see the other email.
> >>>>>>>
> >>>>>>> If we go for 1, (as you have already done ;) ), we cannot create a
> >>>>>>> restricted mixed space using the addition syntax.
> >>>>>>>
> >>>>>>>   mixed = V1 + V2
> >>>>>>
> >>>>>> Do you mean to then restrict mixed? I'm already doing
> >>>>>
> >>>>> Precisely.
> >>>>>
> >>>>> After a second thought I think that doing it in the constructor is
> >>>>> more clean. Makes life simpler. Is this also possible i UFL, I mean:
> >>>>>
> >>>>>   e = FiniteElement("CG",cell,0,"facet")
> >>>>>
> >>>>> Should it?
> >>>>
> >>>> Is that a PyDOLFIN finite element? If it is, then yes, that's how it
> >>>> should work.
> >>>
> >>> I meant UFL.FiniteElement.
> >>
> >> OK. I don't know if that is the syntax. Kristian will know.
> >>
> >>>> Somewhat related, it's a bit annoying in PyDOLFIN that the code is
> >>>> recompiled when the mesh is changed. Can we avoid this? Should one
> >>>> then create a finite element and a dof map with which to initialise
> >>>> the Pythin FunctionSpace?
> >>>
> >>> Do you mean when one has refined a mesh?
> >>
> >> Yes, or just changed from one simulation to the next.
> >
> > Are you sure we get a recompilation, or is the module just loaded from
> > cache?
>
> It's recompiling. It seems to create code for each mesh I use. If I go
> back to a previously used mesh, it's super fast.

Strange, I get this after an instant-clean

   >>> from dolfin import *
   >>> m = UnitCube(10,10,10)
   >>> m2 = UnitCube(5,5,5)
   >>> V = FunctionSpace(m,"CG",2)
   Calling FFC just-in-time (JIT) compiler, this may take some time.
   >>> V2 = FunctionSpace(m2,"CG",2)
   >>>

As can be seen I do not get a recompile. This is all handed by ffc jit and 
JITObject. Particular JITObject.signature.

Put some print statements in there to debug maybe?

Johan


> Garth
>
> >> As far as I understand it is the
> >>
> >>> reason for this that we need to recompute the dofmap.
> >>
> >> In C++ we create a dolfin::FunctionSpace from a dolfin::FiniteElement, a
> >> dolfin::DofMap and a dolfin::Mesh. The FiniteElement and DofMap are
> >> created with a ufc::finite_element and a ufc::dof_map, respectively.
> >>
> >> In the Python interface, in addition to what we have now, we would need
> >> FiniteElement and DofMap objects, with which we could create a
> >> FunctionSpace (together with a Mesh). Then, hopefully we could avoid
> >> recompilation when just the mesh is changed.
> >
> > If a mesh is refined it is the same mesh, and then it would be cool to
> > just reinitialize the FunctionSpace, or the DofMap. Wouldn't this also be
> > a feature in C++?
> >
> >   FunctionSpace::reinitialize()
> >   FunctionSpace::reinitialize(shared_ptr<Mesh> Mesh)
> >
> > or something like that. The first one when the same mesh is being used,
> > for example after a refinement, and the second one when another mesh is
> > being used. Then we do not have to create a new FunctionSpace.
> >
> > Johan
> >
> >> Garth
> >>
> >>> How is this done in the C++ interface.
> >>>
> >>> Johan
> >>>
> >>>> Garth
> >>>>
> >>>>> Johan
> >>>>>
> >>>>>>      V0 = FunctionSpace(mesh, "CG", 3, "facet")
> >>>>>>      V1 = FunctionSpace(mesh, "DG", 3)
> >>>>>>
> >>>>>>      mixed = V0 + V1
> >>>>>>
> >>>>>> and it works. This type of restriction (we'll need to sort out our
> >>>>>> terminology) impacts the dofmap, which I why I think it needs to be
> >>>>>> involved in the creation of the function space.
> >>>>>>
> >>>>>> Garth
> >>>>>>
> >>>>>>> If we want a restricted mixed space we need to use the constructor.
> >>>>>>>
> >>>>>>>   mixed = MixedFunctionSpace((V1,V2),restriction="facet")
> >>>>>>>
> >>>>>>> Johan
> >>>>>>>
> >>>>>>>> Garth
> >>>>>>>>
> >>>>>>>>> or
> >>>>>>>>>
> >>>>>>>>>     V_dg = FunctionSpace(mesh, "DG", 3)
> >>>>>>>>>     V_cg = FunctionSpace(mesh, "CG", 3)
> >>>>>>>>>
> >>>>>>>>>     mixed = V_cg["scalar"] + V_dg
> >>>>>>>>>
> >>>>>>>>> Garth
> >>>>>>>>>
> >>>>>>>>>> Kent
> >>>>>>>>>>
> >>>>>>>>>> _______________________________________________
> >>>>>>>>>> DOLFIN-dev mailing list
> >>>>>>>>>> DOLFIN-dev@xxxxxxxxxx
> >>>>>>>>>> http://www.fenics.org/mailman/listinfo/dolfin-dev
> >>>>>>>>>
> >>>>>>>>> -----------------------------------------------------------------
> >>>>>>>>>-- -- -- -
> >>>>>>>>>
> >>>>>>>>> _______________________________________________
> >>>>>>>>> DOLFIN-dev mailing list
> >>>>>>>>> DOLFIN-dev@xxxxxxxxxx
> >>>>>>>>> http://www.fenics.org/mailman/listinfo/dolfin-dev
> >>>>>>>>
> >>>>>>>> _______________________________________________
> >>>>>>>> DOLFIN-dev mailing list
> >>>>>>>> DOLFIN-dev@xxxxxxxxxx
> >>>>>>>> http://www.fenics.org/mailman/listinfo/dolfin-dev