ffc team mailing list archive

[Anders Logg <logg@xxxxxxxxx>]

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?

> >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).

> >>>>>>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).

--
Anders

Attachment: signature.asc
Description: Digital signature