dolfin team mailing list archive

[Anders Logg <logg@xxxxxxxxx>]

On Thu, Aug 21, 2008 at 09:50:07AM +0200, Martin Sandve Alnæs wrote:
> 2008/8/21  <kent-and@xxxxxxxxx>:
> >> 2008/8/20 Anders Logg <logg@xxxxxxxxx>:
> >>> On Wed, Aug 20, 2008 at 08:43:08AM -0600, Ostien, Jakob T wrote:
> >>>> Hi,
> >>>>
> >>>> I need to be able to assemble an active set of integration points.
> >>>> Essentially, to determine the set I loop over cells and then loop again
> >>>> over the integration points in the cell and determine if that
> >>>> integration point is active with some criteria.  Then I'd like to be
> >>>> able to assemble that set.
> >>>>
> >>>> This is a problem because currently the element_tensor does not break
> >>>> down into integration points, and I need to take derivatives, so the
> >>>> QuadratureElement in FFC is also ruled out.  I suppose I could
> >>>> calculate derivatives and then project on the QuadratureElement, but
> >>>> that seems sort of unclean.
> >>>>
> >>>> It seems like the recent discussion about integrating at a point (on
> >>>> the UFC list) might help me out here.
> >>>>
> >>>> Any other thoughts on how I might go about this?
> >>>>
> >>>> Jake
> >>>
> >>> I don't know yet. I'm thinking about how much of this should go into
> >>> the compiler and how much should be left to the user.
> >>>
> >>> We're working on similar things (and I know Garth is too). I think a
> >>> common denominator would at least be to be able to evaluate a form at
> >>> a given point.
> >>>
> >>> I like Kent's suggestion earlier about adding a new UFC class to go
> >>> along with the other three integral classes, something like
> >>>
> >>>  point_integral
> >>>
> >>> with tabulate_tensor taking the coordinates for a point as additional
> >>> input.
> >>>
> >>
> >> I liked your version with adding a new function to each *_integral class
> >> better.
> >> The reason is that each *_integral object will correspond to one
> >> foo*dx(i) or bar*ds(j)
> >> in the form definition, and each *_integral object will be able to
> >> compute both an
> >> integral over some cell and the corresponding integrand.
> >>
> >> Although a separate point_integral might make sense for something I
> >> don't know about,
> >> matching them to corresponding cell_integral, exterior_facet_integral
> >> and interior_facet_integral
> >> objects will become messy.
> >>
> >
> > One may consider point evaluation also as a form, that was my point.
> >
> > There is an advantage with adding a new function to the classes in the sense
> > that the operators and functions make sense in certain places.
> > Eg. jump can be evaluated on points in the interior of a facet whereas it
> > should not be used on the interior of a cell (note that the jump function
> > does not
> > have the same meaning in eg. a vertex). This coincide (of course) with the
> > the standard usage of dx, ds, and dS. Although with dx, ds, dS one knows that
> > eg jump is not evaluated at vertices because of the choosen quadrature
> > points.
> 
> I see.
> 
> > An advantage of having a separate form object and a separate dx-like
> > symbol in
> > the form language, say dp or delta, is that one may then define forms that
> > are not
> > based on the dx, ds, and dS forms. Put in another way, one can eg. define
> > functions
> > in this way;
> >
> > f = u*sin(x)*dp
> >
> > (although dp does not look nice and I don't have any good suggestion).
> >
> > Kent
> 
> Maybe we should add both ways to UFC?
> Seems to me they cover two different concepts,
> evaluating the integrand of an existing form at a point,
> and evaluating a form that is defined only at a point.
> 
> // Corresponding to "f = u*sin(x)*dp":
> class point_integral
> {
>     virtual void tabulate_tensor(...) const = 0;
> };
> 
> // Corresponding to "a = u*v*dx":
> class cell_integral
> {
>     virtual void tabulate_tensor(...) const = 0;
>     virtual void tabulate_integrand(...) const = 0;
> };
> + same for other integrals.

When I think of this again, wouldn't it be the easiest to just add an
extra function to each integral class (not an extra integral class).
Both form compilers can already generate quadrature code, so we will
get these functions for free since they will be simple versions of
what we already have (take away the quadrature loop).

It will also make it possible to do things like

  a = inner(grad(v), grad(u))*dx

  A = assemble(a, mesh1 - mesh2)

to integrate a form over the difference of two domains:

  mesh1 - mesh2 <--> \Omega_1 \setminus \Omega_2

Some of the cells in \Omega_1 will be removed and some of them will be
cut in two. The assembler could find out which cells are cut and use
quadrature for those cells (by calling the new function).

The above example is my motivation for wanting to evaluate forms at
points.

-- 
Anders

Attachment: signature.asc
Description: Digital signature