Re: [Ffc] evaluate_integrand
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. -- 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
Re: [Ffc] evaluate_integrand
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? 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? I think it is better not modifying UFC at this point. You can derive a class from ufc::foo_integrals which you can pass information on the geometry of polyhedra to compute the integration points. Mehdi > > -- > Anders > ___ > Mailing list: https://launchpad.net/~ffc > Post to : ffc@lists.launchpad.net > Unsubscribe : https://launchpad.net/~ffc > More help : https://help.launchpad.net/ListHelp ___ Mailing list: https://launchpad.net/~ffc Post to : ffc@lists.launchpad.net Unsubscribe : https://launchpad.net/~ffc More help : https://help.launchpad.net/ListHelp
Re: [Ffc] evaluate_integrand
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. > 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. > I think it is better not modifying UFC at this point. You can derive a > class from ufc::foo_integrals which you can pass information on the > geometry of polyhedra to compute the integration points. I don't understand how that helps. Don't you run into the same problem, that the location o
Re: [Ffc] evaluate_integrand
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 ru
Re: [Ffc] evaluate_integrand
On 13/04/10 15:45, 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 OK. 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) OK. We have a member function of the object with which we construct forms that determines the intersection 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. OK Garth -- Anders ___ Mailing list: https://launchpad.net/~ffc Post to : ffc@lists.launchpad.net Unsubscribe : https://launchpad.net/~ffc More help : https://help.launchpad.net/ListHelp
Re: [Ffc] evaluate_integrand
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. 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 gener
Re: [Ffc] evaluate_integrand
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. For this reason, we would like to generate code for evaluating the integrand of a form at an arbitrary point within the cell. I would propose adding something like this to each of the *_integral classes in the UFC interface: /// Evaluate integrand at given point virtual void evaluate_integrand(double* A, const double * const * w, const cell& c, const double* coordinates) const = 0; What about virtual void tabulate_tensor(double* A, const double * const * w, const cell& c, uint num_quad_points, const double * const * coordinates, const double * weights) const = 0; instead? This would allow more of the quadrature code optimisations to be employed. Garth I suspect/hope this should be easy to add to FFC since the current quadrature code must do something like this, but in addition iterate over the points and compute the weighted sum. This would let a user handle the quadrature manually in cases where that is necessary. Would it be possible to add? -- Anders ___ Mailing list: https://launchpad.net/~ffc Post to : ffc@lists.launchpad.net Unsubscribe : https://launchpad.net/~ffc More help : https://help.launchpad.net/ListHelp ___ Mailing list: https://launchpad.net/~ffc Post to : ffc@lists.launchpad.net Unsubscribe : https://launchpad.net/~ffc More help : https://help.launchpad.net/ListHelp
Re: [Ffc] evaluate_integrand
On Tue, Apr 13, 2010 at 08:27:11PM +0800, 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. For this reason, we would like to generate code > >for evaluating the integrand of a form at an arbitrary point within > >the cell. > > > >I would propose adding something like this to each of the *_integral > >classes in the UFC interface: > > > > /// Evaluate integrand at given point > > virtual void evaluate_integrand(double* A, > > const double * const * w, > > const cell& c, > > const double* coordinates) const = 0; > > > > What about > > virtual void tabulate_tensor(double* A, > const double * const * w, > const cell& c, > uint num_quad_points, > const double * const * coordinates, > const double * weights) const = 0; > > instead? This would allow more of the quadrature code optimisations > to be employed. I thought about this but wasn't sure it was general enough for our needs. I'll need to think about it. -- 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
Re: [Ffc] evaluate_integrand
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 nic
Re: [Ffc] evaluate_integrand
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. 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. Since to ev
Re: [Ffc] evaluate_integrand
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
Re: [Ffc] evaluate_integrand
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
Re: [Ffc] evaluate_integrand
On 13 April 2010 11: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
