ffc team mailing list archive

[Anders Logg <logg@xxxxxxxxx>]

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.

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

--
Anders

Attachment: signature.asc
Description: Digital signature