Subject: [Ffc] Fwd: Re: function on EnrichedElement

From: Marie Rognes <[email-masked]>
Date: Tue, 15 Jun 2010 15:35:40 +0200

-------- Original Message --------
Subject: 	Re: function on EnrichedElement
Date: 	Tue, 15 Jun 2010 15:27:44 +0200
From: 	Marie Rognes <[email-masked]>
To: 	Mehdi <[email-masked]>
CC: 	[email-masked], Anders Logg <[email-masked]>

On 15. juni 2010 15:12, Mehdi wrote:
> On Tue, 2010-06-15 at 14:26 +0200, Marie Rognes wrote:
>   
>> On 14. juni 2010 19:57, Marie Rognes wrote:
>>     
>>> On 14. juni 2010 19:37, Mehdi wrote:
>>>   
>>>       
>>>> On Mon, 2010-06-14 at 19:18 +0200, Marie Rognes wrote:
>>>>   
>>>>     
>>>>         
>>>>> On 14. juni 2010 18:10, Mehdi wrote:
>>>>>     
>>>>>       
>>>>>           
>>>>>> Hi Marie,
>>>>>>
>>>>>> I have a function defined on Enriched Element. I want to have access to
>>>>>> the subfunctions defined on this space. How ffc/ufl can be extended to
>>>>>> handle this issue? 
>>>>>>
>>>>>>   
>>>>>>       
>>>>>>         
>>>>>>             
>>>>> Could you give me a piece of sample code?
>>>>>     
>>>>>       
>>>>>           
>>>> The input is something like this,
>>>>
>>>> Elem1 = VectorElement("Lagrange", triangle, 2)
>>>> Elem2 = VectorElement("Lagrange", triangle, 1)
>>>>
>>>> Element = Elem1 + Elem2
>>>>
>>>> u = Coefficient(Element)
>>>>
>>>> # I want to have this functionality
>>>> u1, u2 = split(u)
>>>>
>>>> Thank you in advance.
>>>>
>>>>   
>>>>     
>>>>         
>>> Will see what I can do tomorrow.
>>>
>>> Note that functions on enriched spaces can be a bit treacherous since
>>> the basis is not a nodal basis.
>>>
>>>       
>> I've been thinking some about this. Let me try to explain why this is
>> nontrivial.
>>
>> Say, we have two element spaces
>>
>>     V = span ( {v_i} )
>>     Q = span ( {q_j})
>>
>> and create an "enriched space"
>>
>>     W = V + Q
>>
>> so that
>>
>>     W = span ( {v_i, q_j})
>>
>> We also define the degrees of freedom on W by taking the set of
>> all degrees of freedom on V (L_i) and all degrees of freedom on Q (K_j)
>>
>>     ( {L_i}, {K_j} )
>>
>> At this point, we are starting to tweak the ffc element framework a bit,
>> because the basis for W is _not_ a nodal basis with respect to the
>> degrees of freedom: it might be (actually, quite often is) that
>>
>>     L_i (q_j) \not = 0
>>
>> for all i, j and vice versa for K_j (v_i).)
>>
>> A function f in W can be represented as
>>
>>     f = \alpha_i v_i + beta_j q_j
>>
>> However, note that for some degree of freedom L_k,
>>
>>     L_k(f) = \alpha_k + \beta_j L_k(q_j)
>>
>> Hence, the coefficients in the vector do _not_ directly correspond to
>> the degrees of freedom, unless L_k(q_j) = 0 for all k, j. This means
>> that most kinds of function manipulation on enriched elements give you
>> something different than what you might think you get.
>>    
>> There are however exceptions where things are good: Take for instance
>> continuous linears plus bubble:
>>
>>     V = FiniteElement("CG", "triangle", 1)
>>     Q = FiniteElement("B", "triangle", 3)
>>     W = V + Q
>>
>> Since the bubble basis functions are zero on the boundary of each
>> element, with the above notation
>>
>>     L_k (q_j) = 0
>>
>> for all k, j. This means that you can do
>>
>>     f = Function(W)
>>     g = interpolate(f, V)
>>
>> Then g will be what I think you want from
>>
>>     (g, h) = split(f)
>>    
>> In order to get h, you could do
>>
>>     h = f - g
>>
>>
>> It would be possible to extract the vectors {\alpha_i} and {\beta_j}
>> from a Function on an enriched elements, but this requires quite a bit
>> of work, and is mainly DOLFIN related (rather than FFC/UFL).
>>     
> Actually we were handling this inside main file before, but we didn't
> like it. It makes code dirty especially for the nonlinear problems in
> which we need to obtain these functions in each iteration explicitly.
>
>  We think that handling enriched elements inside FFC/UFL is not
> consistent with the rest. This causes problems while working with
> enriched elements.
>   

So, let's take this discussion on a mailing list instead?

--
Marie

---

From: Mehdi <[email-masked]>
Date: Wed, 16 Jun 2010 13:41:05 +0200

On Tue, 2010-06-15 at 15:35 +0200, Marie Rognes wrote:
>                                CC: 
> [email-masked], Anders Logg
> <[email-masked]>
> 
> 

This is not the case for the enrichment space(Q) in the partition of
unity method. They are often defined on the same dofs that we have
defined standard space V.

> >> of work, and is mainly DOLFIN related (rather than FFC/UFL).

To handle this inside DOLFIN, we extract required data form vectors. Why
performing the approach inisde FFC/UFL is difficult? Isn't it just
enough to pick up the components corresponding to a specific space from
enriched space?

Mehdi

> _______________________________________________
> Mailing list: https://launchpad.net/~ffc
> Post to     : [email-masked]
> Unsubscribe : https://launchpad.net/~ffc
> More help   : https://help.launchpad.net/ListHelp

---

From: Marie Rognes <[email-masked]>
Date: Wed, 16 Jun 2010 14:03:33 +0200

> enough to pick up the components corresponding to a specific space from

I don't see how to do this. If you or anyone else know how, go ahead :)

--
Marie

---

From: Mehdi <[email-masked]>
Date: Wed, 16 Jun 2010 15:51:58 +0200

> I don't see how to do this. If you or anyone else know how, go ahead :)

OK. let me give it a try. Maybe this is a hack, but I think we may need
to attach offsets for the sub functions derived from enriched space. Is
there any reason why not this should work?

 Mehdi

---

From: Marie Rognes <[email-masked]>
Date: Wed, 16 Jun 2010 16:20:21 +0200

> there any reason why not this should work?

I probably misunderstood your original request, thinking that you wanted
to extract
(v, w) from {v + w}...

One way of doing what I now _think_ that you want, is to add some
functionality to FFC/UFL to convert an enriched mixed element to a mixed
enriched element and then do split. This is feasible. But I still think
that this should _not_ be the default behaviour for enriched mixed elements.

--
Marie

---

From: "Garth N. Wells" <[email-masked]>
Date: Wed, 16 Jun 2010 15:45:27 +0100

On 16/06/10 15:20, Marie Rognes wrote:
>>>>>                               From:
>>>>> Marie Rognes<[email-masked]>
>>>>>                                 To:
>>>>> Mehdi<[email-masked]>
>>>>>                                 CC:

A simple example of what's missing is the Stokes mini element. It's not 
presently possible to prescribe only one of the displacement components 
(e.g. u_x) on the boundary. Using V.sub(0).sub(0) in DOLFIN to get the 
u_x component breaks.

Garth

---

From: Kristian Oelgaard <[email-masked]>
Date: Wed, 16 Jun 2010 16:55:00 +0200

On 16 June 2010 16:45, Garth N. Wells <[email-masked]> wrote:
> (e.g. u_x) on the boundary. Using V.sub(0).sub(0) in DOLFIN to get the u_x
> component breaks.

Could it work if instead of doing:

P1 = VectorElement("Lagrange", triangle, 1)
B = VectorElement("Bubble", triangle, 3)
Q = FiniteElement("CG", triangle, 1)
Mini =  (P1 + B)*Q

you did:

P1 = FiniteElement("Lagrange", triangle, 1)
B = FiniteElement("Bubble", triangle, 3)
Mini = MixedElement([P1 + B], [P1 + B])*Q

Kristian

---

From: "Garth N. Wells" <[email-masked]>
Date: Wed, 16 Jun 2010 16:28:45 +0100

On 16/06/10 15:55, Kristian Oelgaard wrote:
> On 16 June 2010 16:45, Garth N. Wells<[email-masked]>  wrote:
>>
>>

Yes, that will work, but is it what we want? It's a bit clumsy and we 
lose the dimension-independence since in 3D we would need

   P1 = FiniteElement("Lagrange", tetrahedron, 1)
   B = FiniteElement("Bubble", tetrahedron, 3)
   Mini = MixedElement([P1 + B], [P1 + B], [P1 + B])*Q

Garth

---

From: Kristian Oelgaard <[email-masked]>
Date: Wed, 16 Jun 2010 18:04:13 +0200

On 16 June 2010 17:28, Garth N. Wells <[email-masked]> wrote:
>>>>>>>>>> At this point, we are starting to tweak the ffc element framework
>>>>>>>>>> a
>>>>>>>>>> bit,
>>> u_x
> Yes, that will work, but is it what we want? It's a bit clumsy and we lose
> the dimension-independence since in 3D we would need
>

Mehdi and I discussed this a bit, one way to get around this in FFC is to let
VectorElement accept a FiniteElement as argument, then you can do

element = VectorElement(V + Q)

and still be dimension independent.

Or in UFL we can tweak the '+' operator, such that enriching a
VectorElement means enriching each of the components of 'self' with
the components of 'other'. For this to work the dimension of the two
vector elements must of course be identical but I guess that will
always be the case, otherwise we throw an error.

Kristian

---

From: Mehdi <[email-masked]>
Date: Wed, 16 Jun 2010 22:49:19 +0200

On Wed, 2010-06-16 at 16:20 +0200, Marie Rognes wrote:
> >>>> Marie Rognes <[email-masked]>
> >>>>                                To: 
> >>>> Mehdi <[email-masked]>
> >>>>                                CC: 

I also have this request. As you said I want to have access to
sub-functions defined on enriched element, i.e. (v, w) from {v + w}
where v and w are scaler.

The other request, which we have discussed more, is handling enriched
mixed element, such that we can have access to the sub elements(Note
that the sub element should be itself an enriched element). 

P = V + Q

where V and Q are vector elements 

Mehdi

---

From: Mehdi <[email-masked]>
Date: Thu, 17 Jun 2010 15:44:23 +0200

> always be the case, otherwise we throw an error.

I will go for this option. This allows us to have simpler code and
preserves accessing to the sub-elements of enriched mixed element.

Mehdi

---

From: Marie Rognes <[email-masked]>
Date: Fri, 18 Jun 2010 01:44:10 +0200

> preserves accessing to the sub-elements of enriched mixed element.

How do you plan on handling elements such as the following (relevant in
connection with the PEERS element for linear elasticity) with this approach?

    V = FiniteElement("RT", "triangle", 1)
    Q = VectorElement("B", "triangle", 3)
    W = V + Q

    v = TestFunction(W)
    u = TrialFunction(W)
    a = dot(v, u)*dx

--
Marie

---

From: Kristian Oelgaard <[email-masked]>
Date: Fri, 18 Jun 2010 11:38:16 +0200

> Q = VectorElement("B", "triangle", 3)

I was planning on throwing an error :)
I did not know that one would ever want to enrich a scalar element
with a vector bubble function, but since this appears to be the case
the only two options, as I see it, are:

1) do nothing such that one will have to do

manually.

2) have VectorElement accept FiniteElementBase as argument to allow

  Mini = VectorElement(P1 + B)*Q

which will reduce the labour of defining elements and still be
independent of dimensions.

We could of course try to 'guess' what the user wants but I don't
think that is a good idea.

Kristian

---

From: Marie Rognes <[email-masked]>
Date: Fri, 18 Jun 2010 12:05:13 +0200

> I was planning on throwing an error :)

Please don't :)

> with a vector bubble function, 

Since RT is a vector-valued element, enriching it with a vector bubble
function
is well-defined.

Strictly speaking, my example is not the enrichment of one UFL VectorElement
with another VectorElement. So, I guess you could overload + for
VectorElement (and TensorElement) only. However, that would make

    V = FiniteElement("CG", "triangle", 1)
    V = V*V
    B = VectorElement("B", "triangle", 3)
    W = V + B

and

    V2 = VectorElement("CG", "triangle", 1)
    W = V2 + B   

behave differently, which I imagine could be rather confusing. 

> but since this appears to be the case

Option 2) sounds neat and useful in other contexts too :)

> think that is a good idea.

Agree.

--
Marie

---

From: Kristian Oelgaard <[email-masked]>
Date: Fri, 18 Jun 2010 12:23:43 +0200

> Since RT is a vector-valued element, enriching it with a vector bubble

Ha, I completely missed the 'RT', that's what you get for working with
'CG' only :)
Now it makes a lot more sense to me.

> VectorElement (and TensorElement) only. However, that would make

Yes, that's what we had in mind, then instead of an error we just
return an EnrichedElement, then it's up to the user to make sure that
the enrichment makes sense. I haven't looked at the code in detail
now, but maybe there are other things we need to check for.

> behave differently, which I imagine could be rather confusing.

They do?
I would expect both instances to end up meaning
W = MixedElement([FiniteElement('CG', 'triangle', 1) + B]*2)

or can you spell it out for me?

> Option 2) sounds neat and useful in other contexts too :)

We of course still do this even if we decide on option 1).

Kristian

---

From: Marie Rognes <[email-masked]>
Date: Fri, 18 Jun 2010 12:34:54 +0200

> Now it makes a lot more sense to me.

Ok :)

>> behave differently, which I imagine could be rather confusing.

Assumption A: If you only overload + for VectorElement and TensorElement.

Under A, yes.

> or can you spell it out for me?

Again, under A:

In the first case, V becomes a MixedElement, and then MixedElement +
VectorElement -> EnrichedElement

In the second case V is a VectorElement and then VectorElement +
VectorElement -> VectorElement

> We of course still do this even if we decide on option 1).

Sorry, this statement was too deep for me: we of course do something
even if we do nothing?

--
Marie

---

From: Mehdi <[email-masked]>
Date: Fri, 18 Jun 2010 12:35:07 +0200

> We of course still do this even if we decide on option 1).

Option 2 is also fine with me. Let's see what Garth thinks.

Mehdi

---

From: Mehdi <[email-masked]>
Date: Fri, 18 Jun 2010 13:20:17 +0200

> Assumption A: If you only overload + for VectorElement and TensorElement.

I think if we want to overload +, by default we should treat these two
cases equally. 

Another advantages of using this approach is we don't need to define
unnecessary scaler elements to just get mixed enriched element. The
approach of defining scaler elements can be annoying, if we want to use
both enriched and non-enriched vector elements(which is the case often
for me). 

If we overload +, it is just enough to have:

V = VectorElement("CG", "triangle", 1)
B = VectorElement("B", "triangle", 3)   
M = V + B

We can use M and V both to define our functions.

If we want to extend VectorElement, this would be,

V1 = FiniteElement("CG", "triangle", 1)
B1 = FiniteElement("B", "triangle", 3)   
M = VectorElemnet(V1 + B1)

V = VectorElement("CG", "triangle", 1)

So we also need to define V explicitly.

Mehdi

---

From: Kristian Oelgaard <[email-masked]>
Date: Fri, 18 Jun 2010 13:43:16 +0200

> I think if we want to overload +, by default we should treat these two

Well, I never intended to overload '+' for just Vector and Tensor elements,
I would overload the MixedElement class which is the base class for
the two special types.
Then the two cases will result in the same element, we should of
course check for equal length of the two mixed elements, and in case
of nested mixed elements, we let the '+' operator handle this on the
sub elements.

> We can use M and V both to define our functions.

This is definitely and advantage.

> M = VectorElemnet(V1 + B1)

But we can still extend the vector element to enable this it you think
it will be useful, that was the essence of my earlier and confusing
remark which should have been:

'We CAN of course still do this even if we decide on option 1)'

Kristian

---

From: Marie Rognes <[email-masked]>
Date: Fri, 18 Jun 2010 14:12:15 +0200

> of nested mixed elements, we let the '+' operator handle this on the

Ok, then what will be the effect for a mixed/vector version of my
previous example? (which is even more relevant for the PEERS element ;) )

    V = FiniteElement("RT", "triangle", 1)
    V = V * V
    Q = VectorElement("B", "triangle", 3, 4)

    W = V + Q

The above W will then be an EnrichedElement of vector-valued elements?

    V = V * V
    Q = Q * Q

    W = V + Q

While this will be a MixedElement of vector-valued Enriched elements?

My main interest is just to keep the + operator "predictable".

> 'We CAN of course still do this even if we decide on option 1)'

Ok, sentence makes sense now!

--
Marie

---

From: Kristian Oelgaard <[email-masked]>
Date: Fri, 18 Jun 2010 15:08:24 +0200

On 18 June 2010 14:12, Marie Rognes <[email-masked]> wrote:
> Ok, then what will be the effect for a mixed/vector version of my previous
> example? (which is even more relevant for the PEERS element ;) )
>

This will (and should) crash for sure, len(V) = 2, len(Q) = 4; so it
makes no sense to enrich the sub elements.

> Q = VectorElement("B", "triangle", 3)

V0 = FiniteElement("RT", "triangle", 1)
V = V0 * V0
Q0 = VectorElement("B", "triangle", 3)
Q = Q0 * Q0
W = V + Q

This is will end up as:

E0 = V0 + Q0 # EnrichedElement([V0, Q0])
W = MixedElement( [E0,  E0] )

and the example we had before:

V0 = FiniteElement("CG", "triangle", 1)
V = V0 * V0
B0 = FiniteElement("B", "triangle", 3)
B  = B0 * B0
W = V + B

will be:
E0 = V0 + B0
W = MixedElement([E0, E0])

> My main interest is just to keep the + operator "predictable".

That looks pretty predictable to me, but there might be other elements
that I'm unaware of for which the logic breaks.

Kristian

---

From: Marie Rognes <[email-masked]>
Date: Fri, 18 Jun 2010 15:19:38 +0200

> makes no sense to enrich the sub elements.

This will not (and I think should not) crash at the moment: both V and Q
has rank dimension 4, and so enrichment makes sense.

However, if you try to progagate the enrichment to sub-elements, you
will run into trouble.
Hence, I assume that you will either throw an error (which I would not
like, since the construction is not erranous) or return just the
original EnrichedElement.

Agree.

> that I'm unaware of for which the logic breaks.

In the vector RT/B example above, the first version becomes an
EnrichedElement while the second becomes a MixedElement. I don't find
that transparent.

I would really prefer to have the conversion from enriched to mixed be
explicit. But I'll stop arguing now ;)

--
Marie

---

From: Mehdi <[email-masked]>
Date: Fri, 18 Jun 2010 15:25:09 +0200

> that I'm unaware of for which the logic breaks.

If we want to address this issue, the same concept should be also
considered for ElementRestriction. What would we expect if we enrich a
vector element with a restricted vector element?

 V = VectorElement("CG", "triangle", 1)
 B = ElementRestriction(V, dc)

 W = V + B

Since this restricted element is not vector element anymore.

Yours,
Mehdi

---

From: Kristian Oelgaard <[email-masked]>
Date: Fri, 18 Jun 2010 16:16:19 +0200

> has rank dimension 4, and so enrichment makes sense.

I think it should, because you will make FFC guess how you want to
enrich your element V because the number of sub-elements is not equal.
Using the same logic as in the first 'RT' example:

Which is an EnrichedElement([V, Q])

What you want should be implemented as:

W = V + (Q * Q)

Then FFC will know what to do because it sees two VectorElements with
two sub elements and it will proceed to enrich sub elements such that:

W = MixedElement([ V0 + Q, V0 + Q ])

which is a vector version of the first example.

> EnrichedElement while the second becomes a MixedElement. I don't find

It is if you use MY definition :)

> explicit. But I'll stop arguing now ;)

So you would prefer if we just implemented

W = VectorElement(V + Q) ?

Kristian

---

From: Kristian Oelgaard <[email-masked]>
Date: Fri, 18 Jun 2010 16:23:45 +0200

> Since this restricted element is not vector element anymore.

Can't we just check, inside the '+' operator, if an element is an
ElementRestriction and then if the restricted element is a
VectorElement
BTW, shouldn't it be RestrictedElement (like EnrichedElement) rather
than ElementRestriction?

Another possibility is to propagate the restriction to the subelements
when creating the restricted element and then return a VectorElement
with restricted sub elements. Did we try that and fail at some point?
I forgot.

Kristian

---

From: Marie Rognes <[email-masked]>
Date: Fri, 18 Jun 2010 16:28:27 +0200

> W = VectorElement(V + Q) ?

Yes. Or

    W = MixedElement(EnrichedElement)

--
Marie

---

From: Marie Rognes <[email-masked]>
Date: Fri, 18 Jun 2010 16:32:35 +0200

> BTW, shouldn't it be RestrictedElement (like EnrichedElement) rather

Yes, that would be more consistent.

> with restricted sub elements. Did we try that and fail at some point?

FFC does do this for its internal representation of restricted elements...

--
Marie

---

From: Mehdi <[email-masked]>
Date: Fri, 18 Jun 2010 17:02:45 +0200

> BTW, shouldn't it be RestrictedElement (like EnrichedElement) rather

Yes, it would be better.

> with restricted sub elements. Did we try that and fail at some point?

I think so. We tried to propagate this, but it didn't work for mixed
element as far as remember. Give it a try one more time. If we don't
have this in place, I can't use new functionalities for enriched mixed
element.

Mehdi

---

From: Kristian Oelgaard <[email-masked]>
Date: Fri, 18 Jun 2010 17:09:21 +0200

>    W = MixedElement(EnrichedElement)

Not sure what that should accomplish though.
The whole point of this discussion is how we avoid doing:

V = FiniteElement("CG", "triangle", 1)
B = FiniteElement("B", "triangle", 3)
W = MixedElement([V + B]*2)

This is more work and dimension dependent. Extending FFC such that:

W = VectorElement(V + B)

is easier and dimension independent.

But in either case, as you also pointed out, one will probably have to
define elements like

V1 = VectorElement("CG", "triangle", 1)

anyway.

Therefore we proposed that adding VectorElements should behave like
adding any other vectors i.e., component by component.
But it looks like it's difficult to agree on a syntax which is fool
proof to achieve this.

Kristian

---

From: Marie Rognes <[email-masked]>
Date: Fri, 18 Jun 2010 17:15:26 +0200

> W = MixedElement([V + B]*2)

I meant

    W = MixedElement(V + B)

> But it looks like it's difficult to agree on a syntax which is fool

Yes. But you are the ones who need to deal with this, and the ones
who will be doing the implementation, so ... ;)

--
Marie

---

From: Mehdi <[email-masked]>
Date: Mon, 21 Jun 2010 12:26:08 +0200

> I probably misunderstood your original request, thinking that you wanted

Any new idea, how this can be handled? This functionality is essential
to compute jump over discontinuity surface. 

One possibility can be attaching indices of underlying sub functions. 
If we have u from {v + w} where v and w are scaler elements. If we
define,

u0, u1 = split(u)

We can attach indices to u0 and u1,

u0 = u<0, m>
u1 = u<m, m + n>

where m and n are the space dimension of "v" and "w" respectively. Can
this be implemented inside UFL/FFC?

Mehdi