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
signature.asc
Description: Digital signature
_______________________________________________ Mailing list: https://launchpad.net/~ffc Post to : ffc@lists.launchpad.net Unsubscribe : https://launchpad.net/~ffc More help : https://help.launchpad.net/ListHelp
