ffc team mailing list archive

[Anders Logg <logg@xxxxxxxxx>]

On Tue, Apr 13, 2010 at 09:22:39PM +0800, Garth N. Wells wrote:
>
>
> On 13/04/10 21:13, Anders Logg wrote:
> >On Tue, Apr 13, 2010 at 09:06:35PM +0800, Garth N. Wells wrote:
> >>
> >>
> >>On 13/04/10 21:00, Anders Logg wrote:
> >>>On Tue, Apr 13, 2010 at 08:19:33PM +0800, Garth N. Wells wrote:
> >>>>
> >>>>
> >>>>On 13/04/10 17:59, Anders Logg wrote:
> >>>>>On Tue, Apr 13, 2010 at 11:31:38AM +0200, Mehdi Nikbakht wrote:
> >>>>>>On Tue, 2010-04-13 at 10:43 +0200, Anders Logg wrote:
> >>>>>>>On Tue, Apr 13, 2010 at 10:31:51AM +0200, Mehdi Nikbakht wrote:
> >>>>>>>>On Tue, 2010-04-13 at 09:45 +0200, Anders Logg wrote:
> >>>>>>>>>On Tue, Apr 13, 2010 at 07:17:28AM +0800, Garth N. Wells wrote:
> >>>>>>>>>>
> >>>>>>>>>>
> >>>>>>>>>>On 12/04/10 23:35, Anders Logg wrote:
> >>>>>>>>>>>On Mon, Apr 12, 2010 at 10:20:13PM +0800, Garth N. Wells wrote:
> >>>>>>>>>>>>
> >>>>>>>>>>>>
> >>>>>>>>>>>>On 12/04/10 21:49, Anders Logg wrote:
> >>>>>>>>>>>>>On Mon, Apr 12, 2010 at 09:34:38PM +0800, Garth N. Wells wrote:
> >>>>>>>>>>>>>>
> >>>>>>>>>>>>>>
> >>>>>>>>>>>>>>On 12/04/10 21:29, Anders Logg wrote:
> >>>>>>>>>>>>>>>On Mon, Apr 12, 2010 at 09:21:32PM +0800, Garth N. Wells wrote:
> >>>>>>>>>>>>>>>>
> >>>>>>>>>>>>>>>>
> >>>>>>>>>>>>>>>>On 12/04/10 21:19, Garth N. Wells wrote:
> >>>>>>>>>>>>>>>>>
> >>>>>>>>>>>>>>>>>
> >>>>>>>>>>>>>>>>>On 12/04/10 20:47, Anders Logg wrote:
> >>>>>>>>>>>>>>>>>>We are doing some work where we need to do run-time quadrature over
> >>>>>>>>>>>>>>>>>>arbitrary polyhedra.
> >>>>>>>>>>>>>>>>>
> >>>>>>>>>>>>>>>>>We (Mehdi and I) do this already (using UFC), so I don't see why a new
> >>>>>>>>>>>>>>>>>function is required. Can you explain why evaluate_tensor is not enough?
> >>>>>>>>>>>>>>>>>
> >>>>>>>>>>>>>>>>
> >>>>>>>>>>>>>>>>I meant 'tabulate_tensor'.
> >>>>>>>>>>>>>>>
> >>>>>>>>>>>>>>>Which function do you call for evaluating the integrand?
> >>>>>>>>>>>>>>>
> >>>>>>>>>>>>>>
> >>>>>>>>>>>>>>We evaluate it inside ufc::tabulate_tensor. We construct our forms
> >>>>>>>>>>>>>>with an extra argument, say an object "CutCellIntegrator", which can
> >>>>>>>>>>>>>>provide quadrature schemes which depend on the considered cell.
> >>>>>>>>>>>>>
> >>>>>>>>>>>>>That would require a special purpose code generator.
> >>>>>>>>>>>>
> >>>>>>>>>>>>What's wrong with that? FFC won't (and shouldn't) be able to do
> >>>>>>>>>>>>everything. Just adding a function to UFC won't make FFC do what we
> >>>>>>>>>>>>do now. We reuse FFC (import modules) and add special purpose
> >>>>>>>>>>>>extensions.
> >>>>>>>>>>>
> >>>>>>>>>>>Exactly, it won't make FFC do what we need, but we could *use* FFC to
> >>>>>>>>>>>do what we need (without adding a special-purpose code generator).
> >>>>>>>>>>>
> >>>>>>>>>>>>>Having
> >>>>>>>>>>>>>evaluate_integrand would allow more flexibility for users to implement
> >>>>>>>>>>>>>their own special quadrature scheme.
> >>>>>>>>>>>>>
> >>>>>>>>>>>>
> >>>>>>>>>>>>We make "CutCellIntegrator" an abstract base class, so the user has
> >>>>>>>>>>>>*complete* freedom to define the quadrature scheme and the generated
> >>>>>>>>>>>>code does not depend on the scheme, since the scheme may depend on
> >>>>>>>>>>>>things like how the cell 'cut' is represented.
> >>>>>>>>>>>
> >>>>>>>>>>>Then it sounds to me like that generated code is not at all special,
> >>>>>>>>>>>but instead general purpose and should be added to UFC/FFC.
> >>>>>>>>>>>
> >>>>>>>>>>>And the most general interface would (I think) be an interface for
> >>>>>>>>>>>evaluating the integrand at a given point. We already have the same
> >>>>>>>>>>>for basis functions (evaluate_basis_function) so it is a natural
> >>>>>>>>>>>extension.
> >>>>>>>>>>>
> >>>>>>>>>>
> >>>>>>>>>>I still don't see the need for 'evaluate_integrand' unless you plan
> >>>>>>>>>>to call it directly from the assembler side (i.e. DOLFIN). Is that
> >>>>>>>>>>the case? Perhaps you can give me a concrete example of how you plan
> >>>>>>>>>>to use it.
> >>>>>>>>>
> >>>>>>>>>Yes, that's the plan. In pseudo-code, this is what we want to do:
> >>>>>>>>>
> >>>>>>>>>   for polyhedron in intersection.cut_cells:
> >>>>>>>>>
> >>>>>>>>>     quadrature_rule = QuadratureRule(polyhedron)
> >>>>>>>>>     AK = 0
> >>>>>>>>>
> >>>>>>>>>     for (x, w) in quadrature_rule:
> >>>>>>>>>
> >>>>>>>>>       AK += w * evaluate_integrand(x)
> >>>>>>>>>
> >>>>>>>>>     A += AK
> >>>>>>>>>
> >>>>>>>>>All data representing the geometry, the polyhedra, mappings from those
> >>>>>>>>>polyhedra to the original cells etc is in the Intersection class (in
> >>>>>>>>>the sandbox):
> >>>>>>>>>
> >>>>>>>>>   intersection = Intersection(mesh0, mesh1)
> >>>>>>>>>
> >>>>>>>>>Eventually, we might want to move part of the functionality into
> >>>>>>>>>either DOLFIN or FFC, but having access to and evaluate_integrand
> >>>>>>>>>function makes it possible to experiment (from the C++ side) without
> >>>>>>>>>the need for building complex abstractions at this point.
> >>>>>>>>
> >>>>>>>>As far as I understood you want to compute the integration points for
> >>>>>>>>polyhedrons inside Dolfin and evaluate_integrands will just compute that
> >>>>>>>>integrand in that specific integration points. If it is the case how
> >>>>>>>>would you determine what is the order of quadrature rule that you want
> >>>>>>>>to use?
> >>>>>>>
> >>>>>>>The quadrature rule would be a simple option for the user to
> >>>>>>>set. Currently we only have one general rule implemented which is
> >>>>>>>barycentric quadrature.
> >>>>>>
> >>>>>>Then if you use higher order elements, you need to update your
> >>>>>>quadrature rule manually. I think it would be nice to compute your own
> >>>>>>quadrature rule inside tabualte_tensor by using the standard quadrature
> >>>>>>rule.
> >>>>>
> >>>>>Yes, but the problem is that the computation of the quadrature rule is
> >>>>>nontrivial and it might be better to do it from C++ than as part of
> >>>>>the generated code, especially when the code depends on external
> >>>>>libraries like CGAL. See BarycenterQuadrature.cpp.
> >>>>>
> >>>>
> >>>>What we plan is that the generated code would provide a 1D scheme
> >>>>(or possible just a polynomial order) and the C++ quadrature object
> >>>>could apply it on a sub-triangulations.
> >>>
> >>>I don't understand. How would the 1D scheme connect to a
> >>>triangulation?
> >>>
> >>
> >>Construct a rule for triangles from the 1D scheme (like FIAT does)
> >>and apply it on each sub-triangle.
> >
> >ok. Then I think I understand. So the input to your specialized
> >version of tabulate_tensor is a list of triangles?
> >
>
> No, tabulate_tensor for a given cell first asks if that cell is
> intersected, and if it is intersected it then calls a function which
> supplies the quadrature scheme. An argument to this function is a
> standard quadrature scheme (which is sufficient to integrate the
> functions on a normal cell). The FFC generated code is unaware of
> how the quadrature scheme for intersected cells is computed - that's
> decided on the the code which we build on top of DOLFIN. It's in the
> C++ code which we have written that the sub-triangulation of a cell
> on either side of the intersecting surface is performed.
>
> We're writing this up for the FEniCS book. Some of it is in
>
>   http://dx.doi.org/10.3390/a2031008
>
> >>>Btw, we have considered subtriangulating the polyhedra to get a
> >>>quadrature rule, but decided to just to barycentric quadrature for now
> >>>to save time, both implementation and run time...
> >>>
> >>
> >>We sub-triangulate. I'm hoping that some of the work can be moved to CGAL.
> >
> >I think CGAL can handle this but we decided to skip it for the
> >moment. It should be possible to make it fast, especially since it is
> >completely parallel (no need to match sub triangulations on
> >neighboring triangles).
> >
>
> Not only fast - sub-triangulation gets tricky in 3D.
>
> >>>>>>>>Since to evaluate integrands you need to tabulate basis functions and
> >>>>>>>>their derivatives on arbitrary points. How do you want to tabulate basis
> >>>>>>>>functions and their derivatives inside evaluate_integrands?
> >>>>>>>
> >>>>>>>That's a good point. We would need to evaluate the basis functions and
> >>>>>>>their derivatives at a given arbitrary point which is not known at
> >>>>>>>compile-time.
> >>>>>>
> >>>>>>Yes, you need to use the tabulate_basis* functions implemented by
> >>>>>>Kristian. Then I am not the only one who is using them. Good for
> >>>>>>Kristian. ;)
> >>>>>
> >>>>>ok, good. Then there is no principal problem of generating code for
> >>>>>evaluate_integrand (if it is allowed to call tabaulate_basis).
> >>>>>
> >>>>>I don't see what the problem is of adding evaluate_integrand. It is a
> >>>>>natural extension (we have evaluate_basis already), it would be
> >>>>>"simple" to implement (Kristian can correct me if I'm wrong) and it
> >>>>>makes generated code useful for assembly over cut meshes (without the
> >>>>>need for writing a special-purpose code generator).
> >>>>>
> >>>>
> >>>>It's not a problem - just good to scrutinise heavily additions to
> >>>
> >>>Agree.
> >>>
> >>>>UFC to avoid bloat. I'm happy for it to be added. We should make
> >>>>clear that it's for testing/special cases since the performance for
> >>>>computing element tensors will be poor compared to tabulate_tensor.
> >>>
> >>>ok. So we have the following proposed additions to UFC:
> >>>
> >>>   evaluate_integrand
> >>>   tabulate_tensor(x, w)
> >>>   evaluate_basis (two versions)
> >>>
> >>>   + quite a few other suggestions as blueprints
> >>>   https://blueprints.launchpad.net/ufc
> >>>
> >>>Perhaps it's a good time to think about UFC 2.0 and take care of all
> >>>the blueprints.
> >>
> >>I think that we're settled on evaluate_basis, so this can go into
> >>UFC-dev. The other needs more discussion (and more work), and I
> >>won't have time for a while.
> >
> >Have we decided on the name for the second version of evaluate_basis?
> >
> >Should it be evaluate_basis_reference (as hinted in one of the
> >blueprints).
> >
>
> Seems OK at first glance.
>
> Have you added 'evaluate_integrand' and 'tabulate_tensor(x, w)' as
> Blueprints?

I have now:

  https://blueprints.launchpad.net/ufc/+spec/evaluate-integrand

--
Anders

Attachment: signature.asc
Description: Digital signature