May I suggest that you start by adding just an exterior_edge_integral,
without any averaging etc?
That's a fairly straightforward matter that will take you through ufl, ffc,
ufc, dolfin
and give you a good understanding of the whole chain.
You'll need at least:
* A new measure type exterior_edge/de in ufl
* A new class exterior_edge_integral in ufc.h (in ffc)
* A new assemble_exterior_edges function in the dolfin assembler
but there will definitely be other details that pop up on the way.
I suggest modelling it after exterior_facet_integral, and searching for
exterior_facet in
the source code will give you some good pointers to where to begin. I also
suggest
making small pull requests during the process instead of waiting too long,
so you can
get feedback early. You can f.ex. make a pull request with just the measure
definition in ufl,
just the new class in ufc.h, without making everything work first.
Martin
On 4 September 2015 at 15:07, Martin Sandve Alnæs <marti...@simula.no>
wrote:
> I lost the list, lets keep it in the loop.
>
>
> ---------- Forwarded message ----------
> From: Martin Sandve Alnæs <marti...@simula.no>
> Date: 4 September 2015 at 14:29
> Subject: Re: [FEniCS] [RFC] Implementing edge integrals in 3D
> To: Lizao Li <lixx1...@umn.edu>
>
>
> On 14 August 2015 at 17:29, Lizao Li <lixx1...@umn.edu> wrote:
>
>> Thank you very much. This is very helpful. I have my responses inline
>> below:
>>
>> On Aug 14, 2015 04:46, "Martin Sandve Alnæs" <marti...@simula.no> wrote:
>> >
>> > Dear Larry,
>> > let me give a quick overview of how the measure domain type corresponds
>> to
>> > other components in ffc,ufc,dolfin to explain why using ds is no the
>> right approach.
>> >
>> > The measures in ufl have a domain type that explicitly means
>> > dx = integral over cells
>> > ds = integral over exterior facets (one-sided)
>> > dS = integral over interior facets (between two cells + and -)
>> > + some more specialized ones, please ignore those.
>> >
>>
>> > Look at <ffc>/ufc/ufc.h to see that each domain type corresponds to a
>> class
>> > class <domain_type>_integral
>> > with a variant of tabulate_tensor(...) as well as a function
>> > ufc::form::create_<domain_type>_integral(subdomain_id)
>> >
>> > In the dolfin assembler, there is
>> > one loop over cells that calls cell_integral::tabulate_tensor(...)
>> > one loop over exterior facets that calls
>> exterior_facet_integral::tabulate_tensor(...)
>> > one loop over interior facets that calls
>> interior_facet_integral::tabulate_tensor(...)
>> > and notice that the various tabulate_tensor signatures are tailored to
>> match the expected data.
>> > (Note: it would be _possible_ to generalize the tabulate_tensor
>> signature,
>> > but that is a significant extra task to take on, so please ignore that
>> route).
>> >
>> > What you want is
>> > 1) a new assembler loop in dolfin that iterates over the edges and
>> collects the relevant
>> > sets of cells and facets, similar to the interior facet assembly but a
>> bit more involved
>> > 2) a ufc::edge_integral class with a tabulate_tensor signature capable
>> of taking the
>> > data the assembler needs to pass to generated code
>> > 3) a ufl Measure domain_type string "edge" exactly matching the
>> ufc::edge_integral class
>> > (the 1-1 mapping simplifies ffc)
>> > 4) a new shortname for the measure, e.g. "dE".
>> >
>> > If you really need 3 different variants of edge integrals, make 3
>> names, e.g.:
>> >
>> > de - arbitrary choice of one of the cell values (min cell index, for
>> example)
>> > dE - average over all the adjacent cell values
>> > dJ - sum over the jumps at all facets around the edge in the
>> right-handed direction (which happens to be the one I care about the most)
>> >
>>
>> The third case is similar to integrating jumps where there is a
>> difference between the interior and exterior part. So I guess I should have
>> both dj and dJ.
>>
>> > So that's my summary of the measure/integral/assembler design.
>> >
>> >
>> > However I have a clarifying question for you:
>> > Do you mean e.g.
>> > a) integral over edge of (avg(f) * avg(g))
>> > or
>> > b) integral over edge of (avg(f * g))
>> > and correspondingly (with jump(f) something like (f+ - f-))
>> > A) sum_facets [ integral over facet of jump(f)*jump(g) ]
>> > or
>> > B) sum_facets [ integral over facet of jump(f*g) ]
>> >
>> > In other words: do you mean jumps and averages of the entire integrand,
>> > or of specific functions and test functions within the integrand?
>> > Because those are not the same, and I think you might be trying
>> > to put too much into the measure definition here.
>> >
>> This is a very good question. Thanks for bringing it up. Currently I only
>> need (b) and (B). But I think the more general strategy capable of dealing
>> with all the cases outline above are more favorable. In this case, I will
>> have two measures:
>>
>> de - exterior edge measure
>> dE - interior edge measure
>>
>> then I will need more modifiers for forms like
>> f['any'] - f evaluated in an arbitrary adjacent cell
>> f['avg'] - f averaged in all adjacent cells
>> f['jump'] - the sum over jumps off all adjacent facets
>>
>> This looks closer to the design for the facet integral. Is this approach
>> more appropriate?
>>
>>
> Sorry for the slow response!
>
> I would advice trying to match these new features closely to existing
> concepts in UFL.
> Adding a new operator that is similar to something existing requires much
> less thinking
> than adding a radically new concept. In particular, currently there is no
> such thing in
> UFL as a 'modifier for a form'. But maybe you mean 'modifier for an
> expression'?
>
> Your suggestions f['foo'] are similar to the existing restrictions f('+')
> and f('-'), so lets stick to the () notation.
>
> Some points:
>
> 1) I think f('any') can be dropped, simply by defining that functions
> occuring in an integrand in a *de integral will use an arbitrary cell.
> The generated code will assume just one cell, and the assembler
> will just pick the cell with lowest index or something like that when
> computing *de integrals.
>
> 2) Perhaps instead of f('jump') we can instead reuse f('+') and f('-'),
> simply defining + and - to be 'current' and 'previous' cell
> clockwise or counter-clockwise around the edge?
> There is a major benefit to this: reuse of the algorithms
> for propagating restrictions for interior facet integrals.
>
> Not sure about how to best handle the average over cells.
>
> Martin
>
>
>
>> > Martin
>> >
>> >
>> > On 10 August 2015 at 21:02, Martin Sandve Alnæs <marti...@simula.no>
>> wrote:
>> >>
>> >> Don't use ds or dx. One or more new measures must be added. I'm back
>> from vacation on Friday, will reply more in depth then.
>> >>
>> >> 7. aug. 2015 12.15 skrev "Garth N. Wells" <gn...@cam.ac.uk>:
>> >>>
>> >>>
>> >>> On 6 August 2015 at 07:49, Lizao Li <lixx1...@umn.edu> wrote:
>> >>>>
>> >>>> Dear all,
>> >>>>
>> >>>> I plan to implement edge integral assembly in 3D in FEniCS.
>> >>>
>> >>>
>> >>> Nice. We have an open feature issue on this,
>> https://bitbucket.org/fenics-project/dolfin/issues/106.
>> >>>
>> >>>>
>> >>>> It is good for a number of things, for example assembling the
>> canoniacl projection for 3D Nedelec edge elements. One issue is that there
>> are more than 3 cells intersecting at an edge in 3D.
>> >>>
>> >>>
>> >>> My concern has been whether we can do edge integration efficiently
>> without clever analysis of the form. For example, if an edge integral
>> doesn't need all data from all connect cells (and there might be a lot of
>> connected cells), will an assembler that gets all data be performant?
>> >>>
>> >>>>
>> >>>> At the level of UFL, my design is to add,
>> >>>> ds('m') - arbitrary choice of one of the cell values (min
>> cell index, for example)
>> >>>> ds('avg') - average over all the adjacent cell values
>> >>>> ds('jump') - sum over the jumps at all facets around the edge
>> in the right-handed direction (which happens to be the one I care about the
>> most)
>> >>>>
>> >>>> Suggestions, hints, and pointers to a good starting point in
>> particular are welcome~
>> >>>>
>> >>>
>> >>> I think we need something other than ds. Perhaps we need to be able
>> to pass the topological dimension to dx. Take a look in measure.py from UFL
>> for background.
>> >>>
>> >>> Garth
>> >>>
>> >>>>
>> >>>> Best regards,
>> >>>> Larry
>> >>>> --
>> >>>> Lizao (Larry) Li
>> >>>> Univeristy of Minnesota
>> >>>>
>> >>>> _______________________________________________
>> >>>> fenics mailing list
>> >>>> fenics@fenicsproject.org
>> >>>> http://fenicsproject.org/mailman/listinfo/fenics
>> >>>>
>> >>>
>> >>>
>> >>> _______________________________________________
>> >>> fenics mailing list
>> >>> fenics@fenicsproject.org
>> >>> http://fenicsproject.org/mailman/listinfo/fenics
>> >>>
>> >
>>
>
>
>
_______________________________________________
fenics mailing list
fenics@fenicsproject.org
http://fenicsproject.org/mailman/listinfo/fenics
