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? 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... > >>>>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. Or perhaps we should wait until after DOLFIN 1.0? I have limited resources at the moment so we might want to postpone this. -- 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
